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The  Hellmann-Feynman  forces  for  H2  and  N2,  and  the 
dipole  moments  for  LiH,  BH  and  CH  are  calculated  with  the 
Multiple-Scattering  Xa  wavefunctions.  First,  the  muffin-tin 
charge  density  approximation  is  used  to  compute  these  two 
quantities.  The  results  show  that  the  internuclear  forces 
for  homonuclear  diatomic  molecules  are  always  repulsive,  and 
the  dipole  moments  for  LiH,  BH  and  CH  are  far  from  the  correct 
values.  Thus,  it  is  necessary  to  use  the  non-muffin-tin  charge 
density  in  such  molecular  property  calculations.  The  main 
difficulty  in  such  non-muffin-tin  calculations  is  to  evaluate 
some  three-dimensional  integrals  over  the  intersphere  region. 
To  achieve  this,  we  express  the  integrands  in  Prolate  Spher- 
oidal Coordinates  and  reduce  the  integrals  into  some  one-di- 
mensional integrals.  The  results  show  that  both  the  force  cal- 
culations and  the  dipole  moment  calculations  improve  a  lot 
over  the  use  of  the  muffin-tin  charge  density  approximation. 
For  Hg,  the  Hellmann-Feynman  force  curve  is  obtained  and  com- 
pared with  the  curve  for  the  derivative  of  the  total  energy, 
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thus  checking  whether  the  Hellmann-Feynman  Theorem  is  valid 
for  the  MSXa  wavef unctions.  For  N« ,  the  force  components 
contributed  by  each  of  the  electronic  states  are  calculated 
at  the  experimental  equilibrium  separation  and  are  compared 
with  the  Hartree-Fock  results.  For  LiH,  BH  and  CH,  the  dipole 
moments  thus  obtained  are  respectively  7.29  D,  -1.24  D  and 
-1.958  D  at  the  experimental  equilibrium  separations.  Also, 
the  derivative  of  the  dipole  moment  for  LiH  at  the  equili- 
brium separation  is  obtained  to  be  1.27  D/a.u.  We  conclude 
that  the  non-muffin-tin  correction  is  important  in  the  force 
and  dipole  moment  calculations.  Any  further  improvement  has 
to  come  from  the  correction  to  the  MSXa  wavefunctions . 


CHAPTER  I 
INTRODUCTION 

In  the  investigation  of  the  electronic  structure 
and  properties  of  atoms,  molecules  and  solids,  we  often 
encounter  the  problem  of  treating  the  exchange-correlation 
effect  properly.  The  Xa  theory  suggested  by  Slater  provides 
us  a  good  way  to  handle  such  a  problem  (see  Slater( 1972a) ) . 
Together  with  the  muffin-tin  approximation,  the  Xa  theory 
can  be  used  to  investigate  the  molecular  system  in  a  prac- 
tical and  accurate  way.  Such  a  method,  called  the  Multiple- 
Scattering  Xa  method  or  the  MSXa  method,  was  first  suggested 
by  Slater(1965)  and  then  was  implemented  by  Johnson( 1966) . 

Many  of  the  results  obtained  with  the  MSXa  method 
show  that  it  gives  better  results  in  the  description  of 
the  one-electron  features  than  do  the  ab  initio  SCF-LCAO 
methods,  and  is  significantly  much  faster  in  calculation 
(see  Connolly  and  Johnson (1973)).  However,  in  many  cases, 
we  find  that  the  total  energy  as  a  function  of  molecular 
conformation  is  inaccurate.  The  energy  curves  of  many  diatom- 
ic molecules,  e.g.,  show  that  the  systems  are  unbound.  It 
was  pointed  out  by  Connolly  (see  Connolly  and  Johnson( 1973) ) 


that  such  a  discrepancy  may  be  mainly  due  to  the  muffin-tin 
approximation,  not  the  Xa  approximation,  in  the  calculation 
of  the  total  energy.  The  current  version  for  calculating 
the  total  energy  assumes  the  muffin-tin  approximation  for 
the  charge  density.  The  result  may  be  improved  by  taking 
the  charge  density,  instead  of  its  muffin-tin  part,  gererated 
from  the  converged  wavefunctions  in  the  calculation.  This 
was  verified  by  Danese  on  the  calculations  of  carbon  and 
neon  molecules  (see  Danese(1973) ) . 

If  the  discrepancy  in  the  total  energy  calculation 
is  indeed  due  to  the  use  of  the  muffin-tin  part  of  the 
charge  density  and  not  the  Xa  approximation,  and  the  result 
does  improve  by  using  the  complete  form  of  the  charge  densi- 
ty, then  it  will  mean  that  the  converged  MSXa  wavefunctions 
are  not  too  far  from  the  true  Xa  wavefunctions  and  it  may 
be  useful  to  employ  the  MSXa  wavefunctions  in  other  molecular 
property  calculations.  In  this  dissertation,  we  shall  inves- 
tigate the  quality  of  these  MSXa  wavefunctions  directly  in 
the  dipole  moment  and  force  calculations  of  some  diatomic 
molecules. 

It  was  proved  by  Slater  (see  Slater(1972b) )  that  the 
Hellmann-Feynman  Theorem  is  satisfied  by  the  exact  Xa  wave- 
functions.  That  is,  the  quantum  mechanical  force  acting  on 
a  nuclear  coordinate  is  equal  to  the  electrostatic  force 
exerted  by  the  other  nuclei  and  the  electronic  charge  cloud 
which  is  formed  by  the  exact  Xa  wavefunctions.  Thus,  it  is 


useful  to  check  the  quality  of  the  MSXa  wavefunctions  by 
testing  if  the  theorem  is  also  valid  for  these  approximated 
Xa  wavefunctions.  Also,  as  it  is  well  known  that  the  dipole 
moment  is  more  sensitive  to  the  wavef unction  than  is  the 
energy,  it  is  useful  to  check  the  quality  of  the  MSXa  wave- 
functions  by  doing  some  dipole  moment  calculations. 

The  difficulty  in  all  these  calculations  is  that 
we  have  to  evaluate  some  three-dimensional  multi-center 
integrals  over  an  awkward  region.  In  this  dissertation,  we 
solve  this  difficulty  by  a  computational  technique  which 
reduces  the  integrals  into  one-dimensional  integrals  whose 
integrands  involve  sine  and  cosine  integrals.  The  main 
point  in  this  method  is  to  employ  the  Prolate  Spheroidal 
Coordinates.  It  is  hopeful  that  this  method  can  be  used  in 
some  other  molecular  property  calculations.  We  shall  dis- 
cuss it  in  detail  in  Chapter  IV  and  Chapter  V. 

In  Chapter  II,  we  shall  give  a  brief  discussion  of 
the  exact  Xa  method  and  the  approximate  MSXa  method.  In 
Chapter  III,  we  shall  discuss  the  effect  of  assuming  the 
muffin-tin  charge  density  in  the  calculation  of  total  energy, 
force  and  dipole  moment,  thus  showing  the  necessity  of  using 
the  non-muffin  charge  density  in  the  calculations.  Chapter 
IV  will  be  devoted  on  the  force  calculations  and  Chapter  V 
will  be  devoted  on  the  dipole  moment  calculations.  Finally, 
we  shall  conclude  in  Chapter  VI. 


CHAPTER  II 

EXACT  Xa  METHOD  AND  MUFFIN-TIN  APPROXIMATE  Xa  METHOD 

2.1  Introduction 
The  discussions  in  this  chapter  shall  mainly  concern 
the  basic  Xa  theory  and  also  its  muffin-tin  approximate  form 
in  the  application  to  molecular  systems.  First,  we  shall 
briefly  discuss  the  basic. Xa  theory  and  its  exact  one- 
electron  Schrodinger  equations,  in  which  no  specific  restrict- 
ion is  placed  on  the  form  of  the  charge  density  and  the 
potential.  Secondly,  we  shall  discuss  the  MSXa  method,  in 
which  the  charge  density  and  the  potential  are  muffin-tin 
averaged.  The  approximate  Xa  one-electron  Schrodinger  equa- 
tions and  their  solutions  are  discussed.  The  quality  of 
the  solutions  of  the  approximate  Xa  one-electron  equations 
is  one  of  our  main  concerns. in  this  dissertation. 

2.2  Exact  Xa  One-Electron  Schrodinger  Equations 

The  Xa  method  was  suggested  by  Slater  and  has  been 
discussed  in  detail  by  many  others  (see  Slater( 1972a) , 
Slater  and  Johnson(1972) ) . 

The  Xa  one-electron  Schrodinger  equations  for  the 


spin-orbitals  u^(r)  can  be  derived  by  the  variational  princi- 
ple from  an  expression  of  spin  up  and  spin  down  charge  den- 
sity: 

p+(?)  =  X  niu*(?)ui(?) 


i  ♦ 


P+(^>  =  X  niui(r)ui(r) 

Li 

where  n^  are  the  occupation  numbers  and  the  summation  is  over 
the  spin  up  and  spin  down  orbitals  respectively.  Also  we  de- 
fine the  total  charge  density  as 

p(r)  =  p+(r)  +  p+(r) 
Then  the  Xa  total  energy  is 

<Exa>  =  Xni  j  u*(r1)(-7?)ui(r1)d?1 

*   (   V"       -2ZP      ,»    .,,-•      .    1    (f     2o(r1)0(r2)*     * 
+  J   £  ]?^fp(rl)drl  +  2    J]    iIi-JLdr!dr2 

+  I  ]   "t(?l)UXa(?l)dfl  +  I   j  e+(fl)UXa(?l>d?l 

+  |E5|E53  (2.!) 

Here  the  nuclear  coordinates  are  Rp  and  their  atomic 
numbers  are  Zp.  Rydbergs  are  used  as  units  of  energy,  atomic 
units  as  units  of  distance.  The  first  term  is  the  kinetic 
energy.  The  second  term  is  the  nuclear-electron  coulomb 
energy.  The  third  term  is  the  electron-electron  coulomb 
interaction  term.  The  fourth  and  the  fifth  are  the  "exchange- 
correlation"  energy  terms  of  either  spin.  The  last  one  is 


the  nuclear-nuclear  interaction  term. 

The  above  expression  is  exact  except  for  the  exchange- 
correlation  term.  .Here  we  define  it  as: 

Uxl,  (?)  =  -9a((3/4Tr)p  +  (?))1/3  (2.2) 

and  similarly  for  U    (r). 
.        xa 

The  Xa  total  energy  is  a  function. of  the  occupation 
numbers  n.  and  a  functional  of  the  spin-orbitals  u..  It  can 
be  shown  that  if  we  vary  the  u^   to  minimize  the  total  energy 
with  fixed  occupation  numbers  n^,    we  can  get  the  Xa  one- 
electron  Schrodinger  equations  as: 

{~V1  +  Vc(p'?l}  +  Vxa(p,?l)}ui(?l)  =  £i  Ui(?l)    (2\3) 

where  A   -2Zn       f  p(^o^ 

Vc(p,?1)  =  h?  §  :  4-  2)-r—   d?  (2.4) 

c    -1   ■  p.i  |r  -R  I     J     r12 


is  the  electrostatic  potential  at  position  r 


and 


r  -6a((3/4TT)pf(r1))1/3    for  spin  up  u± 


V   (P,r-,)  = 
xa   '  1 


->    1/3 
L  -6a((3/47r)p+(r1))      for  spin  down  u± 

(2.5) 

is  the  exchange  correlation  potential  at  r^. 

Since  the  potential  in  the  Xa  one-electron  equations 

is  a  functional  of  the  spin-orbitals  u^,  the  Xa  method  is  a 

self-consistent  method.  We  start  by  assuming  a  potential 

v  +V    Solve  the  differential  equations  and  find  the  spin- 
c   xa " 


Assuming 
Molecule 
Potential 
V  +V_ 


Solve  Eq.(2.3) 
Obtain  u. 


Form  p 


after 


convergence 


Calculate 

<E   (p)>  etc, 
xav  ' 


Generate 

Potential 

V'+V 
c   ex 


Convergence 


yes 


(  Stop  ) 


Fig.  2.1  Flow  Chart  of  the  Conceptual  Exact  Xa  Method 
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orbitals.  From  these  spin-orbitals  we  form  the  charge  densi- 
ty. Then,  by  equations  (2.4)  and  (2.5),  we  calculate  the 
potential  V^+VX{^,  ; which  is  then  compared  with  the  starting 
potential  V  +V„  .  If  they  are  "almost"  equal  to  each  other, 

C    X(X 

the  calculation  is  said  to  have  converged.  Otherwise  we  start 
with  V'+V  '  and  repeat  the  cycle  again  until  convergence  is 
achieved.  Then  we  obtain  the  spin-orbitals  u-,  and  with 
equation  (2.1)  we  calculate  the  total  energy  <E   >.  The 
conceptual  computational  scheme  is  shown  in  figure  (2.1). 

There  are  some  difficulties  in  this  scheme,  namely: 

(1)  There  is  a  parameter  a  in  equation  (2.2)..  In 
atomic  systems,  the  a  value  can  be  chosen  in  several  ways, 
all  of  which  amount  to  almost  the  same  total  energy  (see 
Lindgren  and  Schwarz(1972) ) .  One  of  them  is  to  set  the  a 
value  such  that  the  total  energy  calculated  be  equal  to 

the  exact  Hartree-Fock  total  energy.  Another  is  to  set  the  a 
value  such  that   the  Hartree-Fock  total  energy  calculated 
with  the  Xa  orbitals  is  a  minimum.  Still  another  is  to  set 
the  a  value  such  that  the  virial  theorem  is  satisfied. 
However,  for  molecular  systems,  the  space  may  not  be  as 
homogeneous  as  that  of  the  atomic  systems,  and  we  need 
some  other  prescriptions  to  determine  the  a  value. 

(2)  In  general,  the  potential  is  a  very  structured 
function.  For  atomic  systems,  the  potential  is  spherically 
symmetric  and  it  is  possible  to  separate  the  solutions  in 
spherical  coordinates.  However,  for  molecular  systems,  the 


differential   equations  usually  cannot  be  separated  in  any 
coordinates,  unless  some  approximations  are  made  in  the  form 
of  the  potential. • 

(3)  Suppose  the  differential  equations  are  solved 
and  the  spin-orbitals  are  obtained,  we  still  have  to  calcu- 
late the  potential  by  equations  (2.4)  and  (2.5).  and  to 
evaluate  the  total  energy  by  equation  (2.1).  All  of  these 
involve  three-  and  six-  dimensional  integrals  with  compli- 
cated integrands. 

Due  to  all  these  difficulties,  the  MSXa  method  is 
suggested.  With  the  muffin-tin  approximation  and  the  multi- 
ple scattering  formalism,  the  differential  equations  can 
be  solved  easily  and  the  expressions  for  the  total  energy 
and  the  potential  become  much  more  simple.  We  shall  discuss 
these  in  the  next  section. 

2.3  Muffin-Tin  Approximate  Xa  One-Electron  Equations 
and  the  Solutions 

The  Multiple-Scattering  Xa  method  (MSXa)  was  suggest- 
ed by  Slater  in  1965  and  its  implementation  was  begun  by 
Johnson  in  1966.  There  are  three  basic  characteristics  in 
the  MSXa  method:  the  exchange-correlation  approximation,  the 
muffin-tin  approximation  for  the  charge  density  and  the  poten- 
tial, and  the  multiple  scattering  formalism. 

First,  we  divide  the  coordinate  space  of  a  molecule 
into  three  regions:  (I)  the  contiguous  spherical  regions 
surrounding  the  different  nuclei,  (II)  an  intersphere  region 
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which  is  outside  of  the  spheres  of  region(I)  but  inside 
region  (III),  which  is  the  exterior  of  a  sphere  enclosing 
all  the  spheres  in  region  (I).  The  muffin-tin  approximation 
is  such  that  we  approximate  the  potential  and  the  charge 
density  in  region  (I)  and  region  (III)  by  a  spherical  aver- 
age, and  those  in  region  (II)  by  a  volume  average,  namely:  . 
f(r)  =  5^rf(r)dQ     if  re  I  or  III 


f(r)  = 


Vtt  J" 


f(r)dr      if  re  II 
1 1  J  1 1 
where  f  is  the  potential  or  the  charge  density  function, 

VTT  is  the  volume  of  the  intersphere   region,  and  f(r) 

is  the  averaged  quantity. 

With  the  muffin-tin  approximation,  it  can  be  shown 
that  the  difficulties  mentioned  in  the  previous  section 
can  be  solved,  namely:  the  differential  equations  for  the 
spin-orbitals  are  separable,  the  potential  can  be  computed 
easily,  and  the  calculation  of  the  total  energy  reduces  to 
the  evaluation  of  one-dimensional  integrals. 

Starting  with  the  muffin-tin  charge  density  p"(r), 
we  have  an  expression  for  the  approximate  Xa  total  energy: 

Tf^V  =  £  n±   U*(r)(-V2)ui(r)dr  +  UNN  +  UT 
i     ^ 


where  u.  are  the  spin-orbitals  that  minimize  <Exct>, 
UNN  is  the  nuclear-nuclear  interaction  energy, 
UT  is  the  sum  of  the  nuclear-electron  coulomb  energy, 

electron-electron  coulomb  energy  and  the  exchange-correlation 
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energy 


+  |  dr1p(r1)Jdr2_^ 


.) 
"12 


Ca  =  -|a(3/4TT)1/3 

With  this  expression  for  the  total  energy,  it  is 
shown  in  Appendix   A  that  the  necessary  and  sufficient 
condition  for  <E   >  be  a  minimum  is  that  the  u^  satisfy 
the  following  differential  equations: 

{_v2  +  vc(p,r)  +  Vex(r)}  Ui(r)  -  e±l±(r)  (2.6) 

where  V  (p",r)  and  V   (p",r)  are  the  muffin-tin  average  of 

V  (p\r)  and  Vex(p,r)  as  defined  in  equations  (2.4)  and  (2.5). 

We  want  to  emphasize  that  the  muffin-tin  approximate 
one-electron  equations  (2.6)  are  derived  from  the  fact  that 
the  charge  density  is  of  muffin-tin  form.  There  is  only 
one  approximation  (besides  the  Xa  exchange-correlation  approx- 
imation) involved  in  this  derivation,  and  that  is  the  muffin- 
tin  charge  density  approximation. 

To  solve  the  muffin-tin  approximate  Xa  one-electron 
equations : 

{   -V  +  vXa(P'r>>  ui<r)  =  £iui(r) 
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where  we  write  V  (p,r)+V   (r)  as  V   (p,r),  we  may  use  the 

C  6X  XCt 

partial  wave  expansions  of  multiple  scattering  theory  in 
the  different  regions  of  the  molecule.  Within  an  atomic 
sphere  3  (including  the  outer  sphere  which  we  call  the  0 
sphere),  the  potential  is  spherically  symmetric.  Thus,  we 
may  expand  u^    for  state  i  into  a  linear  combination  of 
products  of  a  radial  function  and  a  spherical  harmonic 
(here  we  use  the  real  spherical  harmonics  for  computational 
convenience) : 

where  cj,!?  are  expansion  coefficients  to  be  determined,  and 
the  radial  functions  R£(rg,£i)  have  to  satisfy  the  follow- 
ing differential  equations: 

{-  4  -^-(1-2-4-)  +  l^^-   +  Vx  (P,r  )}  R  (r   e^ 

=  EiRaCrg.Ei). 

In  the  intersphere  region,  the  potential  is  a 
constant,  Vjj,  and  the  differential  equation  for  u±  in  this 
region  is: 

{  -  V2  +  (Vjj  -  ei)}  u±(r)  =  0 

which  is  just  the  equation  of  a  free  electron.  Hence  we  may 
express  u,.  as  a  multicenter  partial  wave  expansion: 
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ui(r)  = 


if   £t  >V„ 


if  El  <V„ 


where  ^  =  ( j  V-j-  x — e ±  I  )  , 

Tij  is  the  ordinary  Neumann  function, 

J£  is  the  ordinary  Bessel  function, 

k«  is  the  modified  spherical  Hankel  function  of  the 
first  kind, 

i„  is  the  modified  spherical  Bessel  function, 
and    AjJ^  are  the  expansion  coefficients  to  be  determined. 
With  the  requirements  that  the  spin-orbitals  ui  and 
their  derivatives  have  to  be  continuous  on  the  boundaries  of 
the  spheres,   we  obtain  a  secular  equation  from  which  the 

eigenvalues  £j_  and  the  expansion  coefficients  {  C^m  }  and 

iS 
{  A^m}  can  be  determined.  Appendix  B   gives  the  details 

of  the  procedure  to  determine  the  eigenvalues  and  the  expan- 
sion coefficients. 

2.4  Relation  Between  the  Exact  Xa  Method 
and  the  Muffin-Tin  Approximate  Xa  Method 

In  this  section,  we  shall  discuss  the  relation 

between  the  exact  Xa  spin-orbitals  ui  and  the  spin-orbitals 

u-  derived  from  the  muffin-tin  approximate  one-electron 

equation  (henceforth  called  MSXa  spin-orbitals).  First  of 


Fig.  2,2     Flow  Chart  of  the  Conceptual  Exact 
Xa  Method  and  the  MbXa  Method. 
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all,  we  shall  look  at  the  differences  between  the  two  self- 
consistent  cycles  for  the  conceptual  exact  Xa  method  and 
the  muffin-tin  approximate  Xa  method  (figure  2.2). 

Both  SCF  cycles  start  with  the  superposition  of 
atomic  potentials.  In  the  exact  Xa  method  we  try  to  solve 
the  Schrodinger  equations  (2.3)  to  get  ^  directly  from  this 
potential  whereas  in  the  MSXa  method  we  solve  the  Schrodinger 
equations  (2.6)  to  get  u.  after  the  potential  being  muffin- 
tin  averaged.  Then  the  charge  density  p  is  formed  from 
either  u.  or  u..  In  the  exact  Xa  method  the  potential  is 
formed  from  p  while  in  the  MSXa  method  the  potential  is 
formed  from  the  muffin-tin  part  of  the  charge  density,  - 
namely  "p.  Then  we  repeat  the  cycle  again  until  the  conver- 
gence is  achieved.  Then  the  total  energy  is  computed  with 
p  or  p  as  the  charge  density  respectively  in  the  exact  Xa 
method  and  the  MSXa  method. 

As  we  have  mentioned  in  the  previous  sections,  the 
basic  difference  between  the  exact  Xa  spin-orbitals  u.^  and 
the  MSXa  spin-orbitals  u.  is  that  they  minimize  two  differ- 


ent total  energy  functionals,  <Exa(p)>  and  <EXa(p")>.  In 
other  words,  we  have  two  different  hamiltonians :  H,  the 
hamiltonian  for  the  exact  Xa  method,  and  H,  the  hamiltonian 
for  the  MSXa  method.  The  ui  and  u±   are  the  eigenfunctions 
of  H  and  H  respectively. 

Since  we  had  employed  the  approximate  hamiltonian 
in  our  calculations,  an  important  question  should  arise: 
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how  is  the  quality  of  the  MSXa  spin-orbtials  ui?  or  how 
good  are  they  compared  to  the  exact  Xa  spin-orbitals  u.? 
The  calculations  of  the  non-muffin-tin  total  energy  by  using 
the  MSXa  spin-orbitals  u±   done  by  Danese  (see  Danese(1973) ) 
have  given  part  of  the  answer  to  this  question  (we  shall 
give  a  more  detailed  discussion  about  this  calculation  in  the 
next  chapter).  In  this  dissertation,  we  want  to  look  at  the 
quality  of  u.^  in  two  more  aspects:  the  Hellmann-Feynman 
force  and  the  dipole  moment  calculations. 

It  was  proved  by  Slater  that  the  exact  Xa  spin- 
orbitals  (with  the  restriction  that  the  a  parameter  be  a 
constant)  satisfy  the  Hellmann-Feynman  Theorem.  That  is, 
the  electrostatic  force  calculated  with  ui  should  be  equal 
to  the  derivative  of  the  quantum  mechanical  expectation 
value  of  the  hamiltonian  H  with  respect  to  the  coordinate  of 
the  nucleus  upon  which  the  force  is  exerted.  Now,  instead 
of  u.,  we  have  ui;  and  instead  of  the  exact  Xa  total  energy 
<EX  (p)>,  we  have  the  non-muffin-tin  total  energy  <Exa(p)>. 
We  can  check  if  the  u±   satisfy  the  Hellmann-Feynman  Theorem 
and  test  their  quality. 

Also,  dipole  moment  calculations  are  performed  on 
some  diatomic  molecules  by  using  ui  as  the  wavefunctions .  As 
we  know,  the  dipole  moment  is  much  more  sensitive  to  the 
wavef unction  than  the  energy  is.  Therefore,  the  dipole  moment 
calculations  should  be  a  good  check  on  the  quality  of  the 
spin-orbitals  too. 
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Since  we  do  have  some  reasonable  results  for  the 
total  energy  by  assuming  the  muffin-tin  charge  density  "p,  we 
wonder  what  would  be  the  result  if  we  assume  the  same  kind 
of  approximation  in  the  force  and  dipole  moment  calculations, 
So,  before  we  discuss  the  force  and  dipole  moment  calcula- 
tions  by  using  the  MSXa  spin-orbitals  u^,  we  shall  first 
investigate  these  two  kinds  of  calculations  by  assuming  the 
muffin-tin  charge  density  p  in  the  next  chapter.  We  shall 
see  that  the  results  are  not  satisfactory  and  hence  it  is 
necessary  to  use  the  p  formed  from  u. ,  instead  of  the  crude 
approximation  p,  in  such  calculations.  The  distinction  be- 
tween these  two  different  approaches  is  also  shown  in  figure 
(2.2). 


CHAPTER  III 

MUFFIN-TIN  CHARGE  DENSITY  APPROXIMATION 

3.1  Introduction 
In  the  last  chapter,  we  discussed  the  exact  Xa 
one-electron  equation.   Instead  of  solving  this  exact  equa- 
tion, we  solve  the  muffin-tin  approximate  Xa  one-electron 
equation  which  is  obtained  by  varying  the  Xa  total  energy 
functional  by  assuming  a  muffin-tin  charge  density  approxi- 
mation. Since  we  had  assumed  this  approximation,  we  would 
like  to  see  how  good  it  is.  One  way  to  check  the  validity 
of  the  muffin-tin  charge  density  approximation  is  to  look 
at  some  of  the  molecular  property  calculations  using  the 
muffin-tin  averaged  charge  density.  In  this  chapter,  we 
shall  discuss  three  such  calculations:  the  total  energy 
calculations,  the  electrostatic  force  calculations  and 
the  dipole  moment  calculations .' With  the  results  obtained, 
we  conclude  that  in  most  cases  it  is  not  satisfactory  to 
use  the  muffin-tin  charge  density  in  the  molecular  calcula- 
tions. Also,  we  shall  discuss  the  non-muffin-tin  total 
energy  calculations  done  by  Danese.  With  the  success  of 
his  calculations,  it  seems  that  we  should  use  the  charge 
density  generated  by  the  spin-orbitals  ui;  rather  than  the/. 
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muffin-tin  charge  density,  in  the  molecular  property  calcu- 
lations. 

3.2  Muffin-Tin  Xct  Total  Energy  Calculations 
As  we  have  shown  in  the  last  chapter,  we  calculate 
the  Xa  total  energy  by  using  the  muffin-tin  charge  density. 
In  solid  state  systems,  especially  the  close-packed  system, 
the  use  of  the  muffin-tin  approximation  is  usually  quite 
successful.  However,  for  molecular  systems,  it  is  not  uni- 
formly good.  Although  the  absolute  values  of  the  total  ener- 
gy are  quite  close  to  the  experimental  results,  we  have  some 
bad  results  when  we  consider  the  dissociation  energy  and 
the  conformation  of  the  systems.  The  calculation  by  Connolly 
and  Sabin(see  Connolly  and  Sabin(1972))  on  the  water  mole- 
cule led  to  a  minimum  total  energy  for  a  linear  molecule 
rather  than  for  the  bent  form  known  to  be  correct  experimen- 
tally. The  work  in  many  diatomic  molecules  indicates  no 
binding  too. 

It  is  noticed  that  the  results  are  worse  in  loose- 
packed  systems  than  in  close-packed  systems.  It  was  suggest- 
ed that  this  is  because  the  intersphere  region  is  larger  in 
the  previous  case  than  in  the  latter  case,  and  the  muffin- 
tin  effect  is  correspondingly  more  severe.  This  was  verified 
by  the  work  done  by  Danese  and  Connolly*  Danese  performed 
the  calculations  for  the  non-muffin-tin  Xa  total  energy  of 
carbon  and  neon  molecule  by  using  the  MSXa  spin-orbitals 
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rather  than  the  muffin-tin  charge  density  (see  Danese 
(1973)).  The  result  is  dramatically  improved  for  C2 .  One 
finds  binding  very  much  closer  to  the  experimental  value 
than  in  Hartree-Fock  calculation,  and  somewhat  closer  than 
the  result  of  the  configuration  interaction  calculation. 
For  Ne?,  one  also  finds  binding. 

With  the  calculations  on  the  total  energy,  one  may 
conclude  that  the  approximation  of  charge  density  has  quite 
a  large  effect  on  the  molecular  properties.  In  the  later 
sections,  we  shall  show  such  muffin-tin  charge  density  effect 
on  the  electrostatic  force  and  the  dipole  moments. 

3.3  Force  Calculations 

In  the  last  section,  we  have  shown  that  the  use  of 
MSXa  spin-orbitals  ^.causes  a  big  improvement  in  the  total 
energy  calculations.  In  this  section,  we  shall  look  at  u^ 
in  a  more  elaborate  way. 

As  we  know,  the  error  in  the  total  energy  calcula- 
tion is  of  second  order  in  the  error  in  the  wavef unction. 
We  would  like  to  look  at  some  molecular  properties  which 
are  more  sensitive  to  the  wavefunction .  The  force  exerted 
on  a  nucleus  by  the  electronic  charge  cloud  and  the  other 
nuclei  is  a  good  test  of  the  quality  of  the  wavefunctions . 
The  Hellmann-Feynman  Theorem,  which  we  shall  discuss  in 
more  detail  in  the  next  chapter,  provides  a  way  to  calculate 
this  force.  It  was  proved  (see  Slater(1972b) )  that  this 
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theorem  holds  for  exact  Xa  spin-orbitals  u.  (with  the  re- 
striction that  the  a  parameter  be  a  constant).  However,  it 
is  well  known  that  for  approximate  wavefunctions  the  theorem 
often  gives  rather  unreliable  results.  Thus,  it  should  be 
a  good  test  of  the  MSXa  spin-orbitals  by  employing  them  in 
a  calculation  of  the  force  exerted  on  a  nucleus  of  a  mole- 
cule by  means  of  the  Hellmann-Feynman  Theorem. 

We  shall  discuss  such  a  calculation  in  Chapter  IV. 
Before  that,  we  may  wonder  what  the  result  would  be  if  we 
use  the  muffin-tin  charge  density  in  such  a  calculation, 
just  as  in  the  total  energy  calculation.  We  shall  show  in 
the  following  that  for  homonuclear  diatomic  molecules,  the 
force  obtained  with  such  a  muffin-tin  approximation  would 
be  always  be  repulsive.  This  cannot  be  physically  possible 
because  that  would  mean  there  is  no  binding  for  such  a 
system. 

In  the  MSXa  calculation  of  a  homonuclear  diatomic 
molecule,  the  space  is  partitioned  into  three  regions:  two 
atomic  spheres  surrounding  the  two  nuclei  A  and  B,  each  with 
nuclear  charge  Z,  an  outer  sphere  enclosing  the  two  atomic 
spheres,  and  the  intersphere  region  II,  as  shown  in  figure 
(3.1). 

By  symmetry,  the  force  on  nucleus  A  is  along  the  z- 
axis  and  is  given  by: 


J 


2Z  f  P(r)  cos6A  d3r  _  gZ  (3  1} 
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Fig.  3.1  Partitioning  of  the  Space  of  a  Homonuclear 
Diatomic  Molecule  in  the  MSXa  Method. 
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where   p  is  the  charge  density  in  atomic  unit, 

0.  is  the  azimuthal  angle  measured  with  nucleus  A  at 

the  center, 

r»  is  the  distance  from  nucleus  A, 

RAT,  is  the  separation  of  the  two  nuclei, 
AB 

FA  is  in  Rydberg/a.u.  units. 
Now,  if  the  muffin-tin  charge  density  approximation 
is  assumed,  the  quantity  p  in  equation  (3.1)  will  be  replaced 
by  "p  and  the  calculation  for  FA  becomes  much  more  simple. 
First,  we  decompose  the  integral  into  four  different  parts, 
one  over  each  region : 

f  p"(r)  cos0A  3 

?2 "  d  r  =  IA  +  IB  +  Iout  +  In 

A 

p  p(r)  cos0A   3 

where   I.  =      2 d  r 

Ja     rA 

f  p(r)  cos0A   o 


•1 

p(f)  cos0A   3 


p(r)  cos0A   3 
I    =       n d  r 

OUt      ■W        -ri 


Itt  =      Z2 d  r 

11     Jll  rA 

We  shall  discuss  each  of  the  above  integrals  in  the  follow- 
ing: 
Integral  Over  Sphere  A:  I  • 

Since  p"(r)  is  spherically  symmetric  inside  sphere  A, 
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we  have : 


r)  cos0A   3 
I ,  =    I   -2 d  r 


=  "J-  ~*T~  *A   ^  JOCOS0A  Sin9A  d°A   Jo  -d*A 


whgre   R.  is  the  radius  of  the  sphere  A  and  is  equal  to  one 

half  of  the  separation  of  the  bond  distance  R. 

Since  the  integral  over  the  coordinate  0.  is  equal 

to  zero,  we  have  I.  =  0. 
A 

Integral  Over  Sphere  B:  Ig 

To  evaluate  the  other  integrals,  we  need  to  expand    "A 

r^ 
in  terms  of  functions  centered  at  nucleus  B  or  at  the  cen- 
ter of  the  outer  sphere.  Such  an  expression  can  be  derived 
from  an  expansion  theorem  for  Neumann  functions  (see  Appen- 
ix  F).  The  expression  is  as  follows: 

r-4Tr  £  I(*+1,0;1,0;*,0)  ~I~Z^-—   r*  Y°(?  )  . 
i»o  RAB 


Y1<V 


<if  rB  <  RAB> 


(3.2) 

y8(?0) 

0 

(3.3) 

-  4u  Y,   IU-l,0;.l,0;i,0)  R^1  yJ_1(Ra-Hb) 

(1£r0>RA0) 

where  Y^1  are  the  real  spherical  harmonics, 

1(1    ,m  .  ;  X,  ,m  ;  I    ,m    )  are  the  Gaunt    coefficients  for 

1    1    z         2    3    1 

real  spherical  harmonics  and  are  defined  by: 


J 


Ymi(r)  Y™2(r)  Y™3(r)  dft       ,  (3.4) 
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and    R.T-,,  R»^  are  the  distances  between  A. and  B  and 
AB'   AO 

between  A  and  0  respectively. 

Inside  sphere  B,  we  have  rg  <  R^g  and  we  can  use 
equation  (3.2)  for  the  evaluation  of  I  _: 

f   p(r)  cos9A 

I  =    \   —2 dJr 

B     JB     rA 


j       p(r)    (4ir/3)*  IlS|A>    d3r 
Je  r? 


=    (4tt/3)*     J       p(rB)(-47r)X  I(£+l,  0;  1,  0;  1 ,  0) 


Yl+1(RA-Rg)      £      o,-  3 

"FET rB  VrB>    d   r 

RAB 


4  V  ^p+i(Ra_Rb) 

=  (4tt/3)*(-4tt)    2.  I(Jt+l,0;l,0;&,0)      *    \+g 

*  =  <»  RAB 


L'^B*    ^    drB      J  Y?( V    d"l 


)Y°(?g)    dfig   =    (47T)*    5£b 

0   ->      •*• 
x  Y1(RA-Rg) 

I      =    (4tt/3)2    (-4tt)    I(1,0;1,0;0,0)   — — *2  D 
B  RAB 

Ra  _ 

B'      B        B 

i 


x    (4tt)*   J(     p(rD)    rD   drT 


Since      1(1,  0;  1,  0;  0,  0)  =  (  1/4tt) 

and  Yi(ra_RB)    =   Y1(°A   =   7r)    =    (3/4tt)*costt   =   -(3/4TT)1 
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f  6  r) 

therefore   IB  =  4"  J,  ^B^   <*B  '   VV 
KAB 

f*8-         2 
where  QB  =    P(rB)4Tr  rB  drB  is  the  electronic  charge 

inside  sphere  B. 

Integral  Over  the  Outer  Sphere:  IQUt 

In  the  outer  sphere,  rQ  >  Rah,  so  we  can  use  equation 

COS0A 

(3.3)  for  the  expansion  of  — n — .   Then  we  have: 

rA 

t  f      ■  ",-\    COS0A   „3 

'out   -       Jeut    ^   ~^TdT 

t  oa 

=    (4tt/3)*  p(r0)4iT    2^  I(£-1,0;1,0;JI,0) 

x  R^"1  Y°fR   R  )    Y^??}  d3r 
x  RA0   Yo-1(RA"R0)  E+T~  d   r 

ro 

=   47r(477/3)* .  X  I(^1,0;1,0;£,0)  ^q1  Y^-1(RA_^0) 


X«l 


R.»t  r0 


p(r0)   2      10- 

FT  r0  dro  W  dpo 


Since    j  Y°(?0)  d^Q  ».(4ir)*  6£Q 


and   2.  f   0   in  the  summation,  therefore    ^ut  =  0 
Integral  Over  the  Intersphere  Region :  I j  j 

In  the  intersphere  region,  the  muffin-tin  charge 
density  is  a  constant,  so  we  can  take  p  out  of  the  integral: 
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II 


p(r)    cosGa      3  - 
2 d  r   =    p 


II 


IJ 


cos© 


A   ^3 


ZT 


d°r 


We   further  decompose   the   integral    into   three  parts: 


coseA    .3r   _ 
—~2 —  d  r  " 


XI     rA 


A1B+H      r. 


coseA    .3 

— 2 —  d   r      - 


COS0 


A   J3 


d  r 


-'A 


COS0A  q 

rr-^   ddr 


Just   as   for  I.    and   !„,    wehave: 


cos0A      o 


d°r  =   0 


'*      rA 


and 


cos0A 

— o —  ddr   =    ( 


aRr 


4,r2   drB)/R2B  =    ttR^/6 


For  the   first   term,    we  have: 


COS0£        o 

— r> —  d  r 


AtBtll      rA 


=     (47T/3)2        { 


Yj(rA)      3r   + 
— a  r  + 


Ii&l  d3r  , 


"W„*IW      rA  \o-VRa8       A 


=    (4 


Tr/3)*{(-4ir)  }_,    I  (  Z+l ,  0  ;  1 ,  0  ;  I,  0).      *      A      ° 

x  1  ro+2  dro    j  Y?<?o>  .dn0 


£+2 
A0 


*  1 


+4tt     2.     K£-1,0;1,0;£,0)    RaV^-I^A"^ 
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g     ro  -TTT-  dr0 

RA0       rQ 


YO(?0)dfi0} 


=  4n.RA0/3  -  2ttRab/3 
Putting  all  of  them  together 


We  can  express  Fry  in  terms  of  the  total  charge  Qjj 
in  the  intersphere  region  and  the  volume  Vjj; 


Vn  =  (47T/3)  (RAB  -  RA  -  RB)  =  ttR^ 


Qii  Qu 

Therefore   !„  -  -±L  v   r^/2  =    "_ 

11  2  RAB 


Comb 


ining  IA,IB,I  f,    Itt  and  the  nucleus  force,  we  have; 


F.  =  2Z 

A 


1 


P(r) 


COS0A   -a       2   9 
9  -  d-'r  _  2ZZ/R^ 


AB 


=2Z  (IA+  h+    W  +  Ilt>  "2Z2/RAB 


=2Z  (0  +  -J^-  .+  0  +  -?0-)  -  2Z2/R2B 
KAB        ^  KAB 


2Z 


R2~  (QB  +  QII/2  "  z) 
RAB 


Since  for  a  homohuclear  diatomic  molecule,  we  have: 


0   =0      -20-0 
'II   wtotal     ^B   wout 
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where  Q    is  the  charge  in  the  outer  sphere. 
out 

Also,  for  a  neutral  homonuclear  diatomic  molecule, 
Q.  ,  ,=2Z.  Thus,  we  finally  have: 

LuLal 

F  =  -2-Zr-  (Q  ■■+  Q      /2  -  Q   -  Q    /2  -Z) 
A    1r2~  ^b   ^total7     ^B    out' 
AB 

2Z 

-  -$jT   (-Qout/.2  +  Z  "  Z) 

-  -  Z-^  (3.5) 

RAB 

Since  Qout  is  always  positive,  we  have  the  result 
that  the  force  exerted  on  the  nucleus  A  is  always  negative, 
or  away  from  the  other  nucleus.  Thus  we  prove  that  the  force 
is  always  repulsive  for  a  homonuclear  diatomic  molecule  if 
we  assume  the  muffin-tin  charge  density. 

Due  to  such  a  result,  we  see  that  it  is  inadequate 
to  employ  the  muffin-tin  charge  density  approximation  in  the 
force  calculations.  As  in  the  total  energy  calculations,  we 
hope  that  improvements  may  be  achieved  if  we  use  the  charge 
density  p  generated  by  the  MSXct  spin-orbit als,  instead  of 
the  muffin-tin  charge  density.  Indeed,  we  have  some  improve- 
ments as  we  shall  discuss  in  Chapter  IV. 

3.4  Dipole  Moment  Calculations 
Another  way  to  check  the  quality  of  the  wavefunctions 
is  to  look  at  the  dipole  moment  calculations.  As  is  well 
known,  the  dipole  moment  is  very  sensitive  to  the  wavefunc- 
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tion  employed  in  the  calculations.  We  shall  use  the  MSXot 
spin-orbitals  for  such  an  investigation.  Before  that,  we 
first  derive  an  expression  for  the  dipole  moment  calculated 
from  the  muffin-tin  charge  density.  We  shall  see  that  such  a 
charge  density  approximation  is  also  inadequate,  and  we  need 
a  more  elaborate  way  to  evaluate  the  dipole  moment. 

Consider  a  neutral  diatomic  molecule  AB.  We  partition 
the  space  as  we  do  in  the  force  calculation.  Since  we  have  a 
heteronuclear  diatomic  molecule,  sphere  A  and  sphere  B  are 
in  general  of  different  size. 

The  dipole  moment  is  defined  as: 

-*■      f  —-*■-*■      3     -*■ 
y  =  -  j  p(r)r  d  r  +  yN 

where   p(r)  is  the  muffin-tin  charge  density, 

yN  is  the  nuclear  contribution  and  is  defined  by: 
yN  =  ZARA  +  ZBRB. 

For  a  neutral  molecule,  it  can  be  shown  (Appendix  D) 
that  the  dipole  moment  is  independent  of  the  choice  of  the 
origin.  Thus,  for  computational  convenience,  we  can  choose 
the  origin  at  the  center  of  the  outer  sphere.  Furthermore,  by 
the  axial  symmetry  of  the  system,  we  note  that  the  only  non- 
vanishing  component  is  the  z-component .  Thus, 


I    -  -y  3 

y  =  -   J  p(r)z0  d  r  +  yN  • 
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To  evaluate  the  integral,  we  decompose  it  into  four 
parts  again: 

j  p(r)z0  d3r 

f   —  -*■     3      f   —  ■+■      3      I   —  -*-     3 

=   1   p(r)zQd  r  +  J   P(r)zQ  d  r  +   J^  p(r)zQ  d  r 

+  Jn  P(?>z0  d3r 


We  shall  consider  them  one  by  one: 


(1) 


[   P(?)zQ 


d  r 


p(r)  (rAcos0A  -  RAQ)  d  r 


f  A  —     3  —  -*■        3 

=  jo  P(rA)rA  drA  \   cos0A  d^A  -  RA0  J   P(r)  d  r 


=  0  -  R  0   =-R  0 
AO^A      AOvA 


(2)     j   p(r)z0d3r 


5  P(r)  (rBcos0B  +  RBQ)  d  r 


=  j   p(rB)r|  drB    cos0B  d^  +  R^  \       p(?)  d3r 


-  0  +  RB0QB  -  RB0Qb 
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(3) 


P(ro)rocos0Q  d  r 


(4) 


f         P(r)zn  d  r 

J  Out  u 

Jout 

=     J7    p(ro>ro  dro    J   cos0o  dQo    =  ° 

1  Lz° 


P(r)zn  d  r 


=    P 


=    PII    ( 


=    PII    { 


d  r  - 


,2°dr"    <. 


[  z0 

fR*B      3 

rodro 

J  o 

(r   cos        +   R      )    d3r    } 
Jb       B  B  BO 


L  zo 


d3r    ) 


cos0o  dfl0   - 


(rAC°S0A  -  RA0>  d  r 


=  PjI  (  0 .  +  RAQVA  -  RBQVB  ) 


where  V.  and  Vg  are  the  volume  of  sphere  A  and  sphere  B 
respectively. 

Combining  all  the  terms,  we  have: 


p(r)z0  d  r  =  -  RAQQA  +  RBQQB  +  Pn(VARA0  -  V^) 


Since   RAQ  =  (RA  +  RB)  -  RA  =  RB 
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Ho   =    (RA   +   V    "   *B   =   RA 


VA  =      4tt   R3/3 


VB   =      4tt    RB/3 


PII    =   QII/VII    =   QIl/    4irRARB(RA+RB> 

r      _    -v  3 

Therefore  p(r)z0  d  r 

=  -  RBQA  +  RAQB  +  4TTRARB(RA+RB)  T"  <RARB  _  RBRA} 


Q 

^11      ■  4tt    ,3  _3. 


-    "    RBQA   +   RAQB   +  QII(RA   -RB>/3 


Also,      yM  =  ZRRRn  -  Z.R,n  =  ZRR,  -Z.R 


N  "  ^BnBO  ~  ^AnAO    ^B^A   "A"B 
Thus,  finally  we  have: 

y  =  -  ■{-  RBQA  +  RAQB  +  Qn(RA-RB)/3}  +  ZRRA  -  Z^ 

=  RB(QA  +  QII/3  "  V  "  RA(QB  +  QII/3  ~  V 

(3.6) 

We  shall  discuss  the  results  computed  with  the  above 
equation  in  Chapter  V.  There  we  shall  find  that  most  of  the 
results  are  not  satisfactory  when  compared  with  the  experi- 
mental values  or  with  the  results  obtained  by  other  methods. 
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The  reason  is  simple:  the  expression  (3.5)  depends  entirely 
on  the  crude  approximation  of  muffin-tin  charge  density.  We 
hope  that  improvements  may  be  achieved  by  using  the  charge 
density  generated  by  the  MSXa  spin-orbitals  instead  of  the 
crude  muffin-tin  density.  We  shall  devote  Chapter  V  to  such 
calculations. 


CHAPTER  IV 

HELLMANN-FEYNMAN  FORCE  CALCULATIONS 

4.1  Introduction 
In  the  last  chapter,  we  saw  the  necessity  of  using 
the  non-muffin-tin  charge  density  generated  from  the  MSXa 
spin-orbitals  u-  in  the  force  and  dipole  moment  calculations. 
To  carry  out  these  calculations,  we  need  to  evaluate  a  three- 
dimensional  integral  whose  integrand  is  a  multi-center  func- 
tion. The  basic  difficulty  of  such  calculations  is  the  eval- 
uation of  the  integral  over  the  intersphere  region,  which  is 
in  general  very  awkward.  In  this  chapter,  we  shall  develop 
an  analytical  way  to  perform  such  integrals  for  the  case  of 
diatomic  molecules,  and  we  shall  apply  it  to  the  force  calcu- 
lations. 

4.2  Prolate  Spheroidal  Coordinates 
Before  we  go  into  the  force  and  dipole  moment  calcu- 
lations, we  first  introduce  the  Prolate  Spheroidal  Coordinates 
that  are  essential  to  our  scheme  of  evaluating  the  integral 
over  the  intersphere  region. 

Let  A  and  B,  with  separation  of  R,  be  the  two  foci 
of  our  coordinates.  Also,  let  (rA,  6A,  4>A)  and  (rB,  9B,  <J>B) 
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B    ^  <(,=(|)A=(j)B 
>ZB 


€ 


Fig.  4.1  Representation  of  a  Point  in  the  Atomic  Coordinates, 
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be  the  two  spherical  coordinates  with  centers  located  at  A 
and  B  respectively.  0,  and  0B  are  both  measured  from  the 
major  axis  in  the  same  direction,  as  shown  in  figure  (4.1). 
The  Prolate  Spheroidal  Coordinates  are  d. fined  as: 


5  =  (rA  +  rB)/R 


n  =  (rA •+  rB)/R  (4.1) 


The  reciprocal  relations  are; 


rA  =  R(C  +  n)/2  rB  =  R(£  -  n)/2 


cosG     =  (l  +  £n)/(?  +  n)         and         cosGg  =  -  (l  -  £n)/U  -  n) 

sineA  ■-  i<g2-  Uil  -n2))'  SineB  -  ^  1)(1  ^»* 

A  5  +  n  B  5  -  n 

(4.2) 

Also,  it  can  be  shown  that  the  volume  element  is: 


d3r  =  (R/2)3(?2  -  n2)  d?  dn  d<j>  (4.3) 


If  the  integral  is  over  the  entire  space,  the  inte- 
gration limits  would  be:  1  S  £  <  °°,  -1  £  n  ^1,  and  0  ^  <J>  $2ir  . 
However,  if  we  are  evaluating  an  integral  over  a  finite  space 
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only,  such  as  the  integral  over  the  intersphere  region,  the 
integration  limits  have  to  be  modified.  We  shall  use  the  two 
atomic  centers  or  one  atomic  center  together  with  the  center 
of  the  outer  sphere  as  the  two  foci,  depending  on  the  form 
of  the  integrands,  in  our  treatment  of  the  integrals. 

4.3  Hellmann-Feynman  Theorem 
It  was  proved  (see  Slater  (1972b))  that  the  exact  Xa 
spin-orbitals,  with  the  restriction  that  the  a  value  be  a 

constant  everywhere,  rigorously  obey  the  Hellmann-Feynman 

9<EXa> 

Theorem.  In  other  words,  the  force  component  -  — — 

9Xp 

acting  on  a  nuclear  coordinate  Xp,  where  <EXa>  is  the  exact 
Xa  total  energy  of  the  system  as  calculated  quantum-mechani- 
cally  from  the  exact  Xa  spin-orbitals,  is  the  force  which 
would  be  calculated  by  pure  classical  electrostatics  acting 
on  this  nuclear  coordinate,  resulting  from  the  electronic 
charge  cloud,  formed  by  the  exact  Xa  spin-orbitals,  and  all 
other  nuclei.  We  can  express  the  theorem  in  the  following 
form: 

_  1<I2^1  -   [  Pxa(?)     *  (    *5e->  d3r  -   £  4  (|2zp;V 
dX  J    xa        ^  l?-RPl  tTP    3Xp   |Rp-Rql 


(4.4) 


where  Z   is  the  atomic  number  of  the  p   nucleus, 

<E„  >  is  the  exact  Xa  total  energy  in  Rydberg  units, 
xa 

p  a  is  the  electronic  charge  density  formed  by  the 
exact  Xa  spin-orbitals, 
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R ■   R  are  the  nuclear  coordinates, 
p'   q 

In  the  above  theorem,  we  want  to  emphasize  that  the 
spin-orbitals  and  the  total  energy  are  those  derived  from 
the  exact  Xa  method.  As  we  have  mentioned  in  the  previous 
chapters,  the  existing  MSXa  calculations  are  not  exact  Xa 
calculations.  They  assume  the  muffin-tin  approximation.  With 
such  an  approximation,  we  can  ask  the  following  two  questions: 

(1)  Does  the  Hellmann-Feynman  Theorem  hold  for  the 
approximate  Xa  spin-orbitals  ui?.If  not,  then  how  severe  is 
the  discrepancy?  In  other  words,  we  want  to  see  how  large 


is  the  difference  between  -  22—   and  Fp(pXa),  where  <Exa: 

3Xp  ^ 

is  the  Xa  total  energy  obtained  with  MSXa  spin-orbitals  ui 


r\, 


ant 


d  Fp(pxa)  is  the  electrostatic  force  exerted  on  the  nuclear 


coordinate  X  by  the  other  nuclei  and  the  electronic  charge 

a. 


a. 


cloud  pxa  formed  from  u^ 

(2)  How  well  does  Fp(pxa)  compare  with  the  experimen- 
tal value  of  force  acting  on  the  nuclear  coordinate  Xp? 

In  order  to  answer  the  above  questions,  we  made  some 
force  calculations  on  homonuclear  diatomic  molecules.  The 
reasons  that  we  pick  the  homonuclear  diatomic  molecules  are 
that  the  a  values  in  these  systems  are  constant  throughout 
the  whole  space  which  is  the  restriction  of  the  Hellmann- 
Feynman  Theorem,  and  that  the  computations  are  easier  and 
the  results  can  be  compared  with  the  experimental  values. 
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4.4  Method  of  Hellmann-Feynman  Force  Calculations 

To  investigate  the  Hellmann-Feynman  Theorem,  we  have 
to  calculate  two  quantities,  the  total  energy  and  the  elec- 
trostatic force  exerted  by  the  electron  cloud  as  expressed 
by  the  integral  on  the  left  hand  side  of  equation  (4.4).  The 
calculation  of  the  total  energy  <EX  >  can  be  done  by  using 
Danese's  program  (see  Danese  (1973)),  as  we  have  discussed 
in  Section  2  of  Chapter  III.  The  main  quantity  that  we  want 
to  evaluate  in  this  discussion  is  the  electrostatic  force 


F  (I      )  = 
e  Kxa 


a     3Xp    |i-Rpi 


Since  p(r)  =  2-  niui(r)ui(r)  where  n^  is  the  occupa- 
tion number  of  the  i   state,  we  can  express  F0(p„  )  as: 


For  a  homonuclear  diatomic  molecule,  the  expression 
for  the  z-component ,  which  is  the  only  nonvanishing  part, 
of  the  force  exerted  on  nucleus  A  is  reduced  to: 

Fe(Pxa)  =  2  ZA  L  ni   «i^Ui(r)V  d3r 
i      v  rA 

To  evaluate  the  above  integrals,  we  partition  the 
region  of  integration  into  four  parts:  one  over  sphere  A, 
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one  over  sphere  B,  one  over  the  outer  sphere,  and  one  over 
the  intersphere  region: 

]  ^(r^C?)  ^  d'r  -  Fj  ♦  4   ♦  F0Jt  ♦  Fj,     (4.5) 

J  rA 


,i    f  %*-*-,.%  ,-*-.  cosOa   ' 
i  =   •  ui(r)u.  (r)  — n— s  d' 

J    JJ  rA 


where  Ft  =  |,  uH(f)u^(f)  r~~"rt  d3r  ,       j  =  A,  B,  out,  II 


We  shall  discuss  each  of  the  integrals  in  the  follow- 
ing: 
Integral  Over  Sphere  A:  F« 

Inside  sphere  A,  we  have 


Ui(r)  =  A  C^  Rfc(rA)£i)  Y£(rA) 


L» 


iA 
where  C   are  the  C-coef f icients, 

7-m/ 


Y'p'Cr.)  are  the  real  spherical  harmonics  centered  at  A, 
Since  we  have  used  real  functions  as  our  basis,  all 
the  C-coef f icients  are  real,  and  ui  =  u . .  Therefore, 


■i     f   a,*  -*■  %      -*■  COS0A  "3 

F^  =     ui(r)ui(r)  _rA  ddr 

JA  rA 

=   Z  Z  CJA   cjA  [  {  R*  (rA)Y»  (rA) 


m  -»-   cosGa  ,   ^ 
(r.)Y0  (r.)  —A  }>dd3 


x  R 


I   I 


^iA   „iA 
c    c 


A 


Ril/rA)^2(^)-?^  drA 


A 
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Y"I(fA>Y™;(?A>COSeA  dt!A 


Z  H  c"micik(4*/3)il(vra.;,l2'"v;1'0> 


*  J.\, 


<  WrA>  drJ 


where  I(l1 ,mi;H2,mz;H    ,m3)   are  the  Gaunt   coefficients  as 
defined  in  equation  (3.4). 

Since  I  (_l1  ,mj  ;  l2  ,m  ;  1,0)  =  0  unless  the  arguments 
satisfy  the  following  conditions: 


£,+'£  +1  is  even 
I   -I    I  <  1  <  A,+A, 

1    *  '  12 


therefore     F*   =         T!  V     C^A      C^A        (4tt/3 


x    i(ai  ,.m;£2,m;l,0) 


R£1(pA)R£se(rA)    drA 


By  the  axial  symmetry  of  the  diatomic  molecule,  the 
C-coef f icients  are  nonzero  only  if  m=mi  for  state  i  (hence- 
forth called  m) .  Therefore,  we  have: 


F.  = 


X 


cjA   cjA    (47T/3)1 


KZ^m-.Z^m;!^)   J   R£  j  (rA)R£z  (rA)dr 


A 
(4.6) 
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Integral   Over   Sphere  B :    F 

. y 

Inside   sphere  B,    we  have 


»!<*>    "      Z    CZ   »i<*B>    *?<*„> 


m    .  ->-    .  _      ,      .  „m 


s°  'b  ■  ,2  <£„  cj»n  j^vvv^vvvv 


cosGa      3 
x    — ? —  d  r 
r^ 
A 

2 

By  equation  (3.2),  we  can  expand  cosOA/r.  in  terms 

of  functions  centered  at  B: 


cosQa         i  Yi(rA) 
__A  =  (.47r/3)*   *  2A 

A  A 

=  -  4tt  (4tt/3)*   2-   J(  M-1,0;1,0;A,0  ) 

£+2     .  B   Jt   B 
.     RAB 

Therefore      F*  =   T  C^C^C-^K^/S)* 

x  T  I(^1,Q;1,0;£,0)  IllliJAJB) 
*s°  RAB 

R£1<rB)R£2(rB)rB+2  drB 


r 


*!,<VY?,<fB>Y?<V.',aB 
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±*K 


K  =       Z      CoB™CoBm(-4Tr>(4Tr/3>*       2]      {I(A+1,0;1,0;A,0) 
B  £-J        lira   lzm  tfnC.,tl 

x      I(£      m;£, ,m;£,0)    ^H^aJbJ 

KAB 

x      j\i(rB)R£2(rB)r^2drB}  (4.7) 


Integral   Over   the  Outer   Sphere:    Fout 

Inside   the  outer   sphere,    we  have: 


Thus      Fout=      S     Ci'?m4°m.jjR£1(r0>YT/?0>R£2(r0>Y?2^0> 

x    }dJr 

rA 

o 
By  equation  (3.3),  we  expand  cos0A/r^  around  the 

center  of  the  outer  sphere  0.  After  integrating  over  n  , 

we  have : 


Ijtl* 


F*   =■  4 
out 


^(4^/3)*   £  cj°racj°m  £  {.l(£'-l,0;l,0.;£-.f0) 

x  iai,m;£2,m;£.',0)  RAQ~  Y^^i^-R^) 

/*R.  (r0)R.  (rQ) 
•/e       rn 
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Integral  Over  the  Intersphere  Region:  Fjj 

In  the  intersphere  region,  u.  can  be  expressed  as: 


r       Y      Ajin   A  j.         A      £j      A7  £      £2m '£2 


u^r)   = 


Z-AjB„n.    («  )^'(J   ) 


f  O^"^",^ 


if  El  >  v„ 


I 


A£   mV(lcrA)Y£    <rA>    +     f  A£   mk£    (^rB)Y£    (rB) 

11  1  /*    '        2  -  2  2 


2> 


10       •  /  NTT113      /*       N 

i      (<rrt)Y      (r    ) 


£,m   £. 


0" 


if   ei   <   Vn 


For  diatomic  molecules,  the  eigenvalues  e±    for  most 
states  (except  the  core  state)  are  greater  than  the  constant 
potential  Vjj  in  the  intersphere  region.  We  shall  look  at 
the  positive  energy  case  (£.>VTI)  in  more  detail  (the  other 
case  can  be  treated  in  a  similar  way).  For  core  states,  the 
electronic  charge  is  concentrated  around  the  atomic  centers 
and  the  contribution  due  to  the  intersphere  region  can  be 
neglected  in  the  core  states.  Thus,  we  have  the  following 
expression  for  Fjj: 

i  1      i/*  -*  ^    ,+      cosOa     ,3 

fii  =  JJt  V"ui(r)  -^  d  r 

-  ^  Ai% am>  'L'1. ■*»•'■> +  £,/£»  AZ*  i'^'1"'1 

+      Z    4?m  *"m  4o<  Wm>'  +   2I.A"m  4!m  IaB«..^.») 
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r  iA   iO   i  V  iB   iO   i 

(4.9) 
where   I^U,  ,£2  ,m)  =  j  ^^(kt^^t^^kt^^t^ 

x   ^Ad3r 
*A 

■I^U,,*,.*)  -  |  ^i(<rB)Y^(?B)n£2(<rB)Y^(rB) 

coseA  ~ 


'2 


IaB(*i, *..">>  "  f  %1<KrA)Y«',<;A)*1«2UrB)Yj'2<";B> 


Iio(1"t"n)  =  (In*1<"A>T?.<VJt.("o)T?,<?o) 


cosGa   ,3 


IB0(£i,£2,m)  =  j  n^(^rB)Y£i(rB)j£2(Kr0)Y^(r0) 


COS0A   ,3 


4 


4(^-l2'm)  =  fu^/Kro)Y^1(?o)^2(Kro)Y^2(?o) 


cosGa^ 
rA 
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To  evaluate  each  of  the  above  integrals,  we  employ 
Prolate   Spheroidal   Coordinates  as  introduced  in  Section  2 
at  the  beginning  of  this  chapter.  First  of  all,  we  want  to 
see  what  the  limits  of  integration  are  for  such  coordinates. 

For  the  intersphere  region,  we  have  the  following 
restrictions  on  £  and  n: 

(1)  rA  +  rB  >,   R  or       5  =  (rA  +  rB)/R  *'  1 

(2)  I  rA  -  rB  |  a  or    |  n  |  =  |(rA  -  rB)/R|  $   1 

(3)  rA  >  RA  =  R/2  or   £  +  n  ^  1 

(4)  rfi  ^  RB  =  R/2  C -  n  >   1 

(5)  From  figure  (4.2),  we  have: 


r2  =  r\   +  (R/2)2  -  2rA(R/2)cos©A 


=  (R/2)2(£+n)2  +  (R/2)2  -  2(R/2)(?+n)(R/2)(l+5Ti)/(C+n) 


(R/2)2  (E2   +  n  -1) 


2     2 
Since   Tq  $  R,  we  have   £   +  n   $.5 


Putting  all  of  the  restrictions  together,  we  have 
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-1 


n=i-S 


Fig.  4.3  Boundary  of  the  Intersphere  Region  in 
Prolate  Spheroidal  Coordinates  with  A 
and  B  as  the  Foci. 
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the  transformed  intersphere  region  as  defined  by:  (figure 
4.3) 


-f(0  $  €  $  f(U 


and 


1  $  £  $  /5 


where 


r  5  -  i 


f(0  = 


if  U<  U  2 


(5  -  5  )f      if  2  ^  C  $  ^ 


Before  discussing  the  evaluation  of  the  integrals, 
we  note  that  all  the  integrands,  except  for  the  factor  of 
the  spherical  harmonics,  are  cf>  -  independent ,  and,  because  of 
the  cylindrical  symmetry  of  the  intersphere  region,  all  the 
integrals  can  be  reduced  into  two-dimensional  integrals 
through  the  (^-integration,  such  as : 


II 


,    ._,.m  ,-*■    .    ,  .  N„m  .-*■    .       2   3 
ri^(KrA)Y£1(rA)n£2(<rB)Y£z(rB)cos0A/rA  d  r 


2.3 


n^  (krA)P£  (cos0A)n£  (xr^P^  (cos0B)cos0A/rAd  r/d<|> 


where  Pm  are  the  normalized  Legendre  polynomials. 


To  evaluate  the  integrals,  we  also  introduce  the  fol- 

m 
lowing  expressions  for  P.,  n.,  and  j,:   • 


l-ifi-w}] 


Pp(cos0.)  =  N„(sin0.)    2-1   %  (cos©.) 


H-m-2v 


A'    *->   -v 
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n»(z)  =  (-1)  +1z_1  {Q°U+i,z)cos(z+£JU)-  Q1  (  £+§,  z)sin(z+£JU) } 


j£(z)  =  z-1  {Q°(£+|,z)sin(z-J^7r)  +  Q1  (£+|  ,z)cos(z-§£tt )} 


m    .   m  (2&+l)U-m)!  » 
where  N,,  =  (-1)   {    2  (£+m) ,  } 


,  Am  _     y     (2£-2v)! 

W      =   .1-1) 

V         2i(£-m-2v)!U-v)!v! 

L4^]     k  -2k 

Q°(*+*,z)  =   2-,  (-l)K(£+i,2k)(2z)   *  (4.10) 

fc=o 


Q^H.z)  =   5]  (-l)k(^+i,2k+l)(2z)" 


<*♦*;»  ■  i^ir 


and    [  b]  =  the  integer  part  of  the  real  number  b. 

With  the  above  expressions  and  some  algebra,  we  have 
the  following  results  for  IAA,  Igg,  IAg,  ^0'  IB0'  and  ^-00 
as  required  by  equation  (4.9): 


2   *3 

d  r 


if  m  -*■  m  -*■         2 

(1)  XAA:     ^1(KPA)yA1(rA)T1£z(lcrA)Y£2(rA)cosVrA 


This  is  a  one-center  integral.  With  some  straight- 
forward algebra,  we  have: 
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{/1v2m(-i)rc,;(KR)r1+^-2m+1-2v+2rc'1+^-2m+1-2v+2r 


,  _  -s   s     2.£i+£2-2m+l-2v+2r-s 
x(kR)    £   (l-£  ) 


x  {4^g£  £2(£i+£2-2m-2v+2r-s+4>S) 


-  (2/kR)  g*  „  (Jl1+£,-2m-2v+2r-s+3,C)}}  (4.11) 


where  v  =  vx+  v2 


c?  =      m! 


'r    (m-r)!r! 


Here  the  g1  -functions  are  defined  as: 


'*i*2   3        <^V  UT.      n*° 


x  A(u'v)  p[U;V'1)(£3+2n+u+v,0 

where  the  A-coef f icients  are  defined  in  terms  of  the  (£+£,k) 
coefficients  through  the  following  relations: 

•JT  A<U'v)(2z)  n  =  QU(llH,z)QV(l2H,z) 

n  =  o 
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and  P^u'v,1^(£, O   are  defined  in  terms  of  the  sine  and 
cosine  integrals  as  below: 


Pixil'1)(£/C)  =  {cos*^i+£2)^  <6u+v.0-fiu+v,2) 

\         cos  t 

+    {sin|(£1+£2)7r    ( 6u+V)2-<Su+v,0) 

i     I                       sin   t    Ju_ 
+   cosi(£   +£,)tt    6   ^      ,}  —  dt 

1         2  U+V,  1  „\        '       t* 

+    {cosiCVV*    («u+vi0+Vv,2) 

,K«(1*#P) 

+   sinKVVir    «u+ViiC«vo-6vl>}    J^^     F 

(2)    :BB:      Jxi^1^rB)Y£i(rB)n£2(KrB)YA2'(rB)coseA/rAd.r 

This   is   a   two-center   integral.    The   result    is   as 
following: 

4 _-«"    <      C-l)t'+^+1(E/2)  «<   L       ill 

BB  £i       J,  2  ^i  ^To       ^j»      *so      t=« 


£,m      i,m     m  £   +Jl2-2m-2v+3r      Jl,+i2-2m-2v+2r 

v;  <  cr  '  ^  n 


xr^  a):2"1  c;  <-ir*- ~2  cs 


x    (KR)S    2(1-C2)S    ^1  +  ^-2m-2V+2r-S    (dH2)    gJi£2(s+l,C) 
-    (5/kR)    g2    £„(s,0}>  <4-12> 


jX-2 
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where  v  =  vi  +  v2 

and  the  g2 -functions  are  defined  as: 


x  A(U'V)  P(u'v'2)(s+2n+u+v+lfO 
n      SL1Z2 


where  the  A-coef f ieients  are  defined  as  in  IAA>  and 
,p!u!V|2V  £)  are  defined  in  terms  of  the  sine  and  cosine 
integrals  as  below: 

kRQ-W) 


cos  t 


r 

+  sinK^+M**     }.        ~ 2~T  dt 

U+V'1-  Vj.^))(2K«-t)2t£ 


+  Uini(l1+£2)ir    (5u+Vj2-5u+v0) 


sin   t 


f 

+  cosi(£  1+^)^6   A      -}  — —   2   A   dt 


+  {cosi(£i-£2)TT    (fiu+y^-^+v^^ 
+   sin£(*   -£    )tt    (6    _-5    .)6   A      .. } 

2V      1         2  '  yO        vl        U+V, 1 


I-dt 


<3>  4:  J„  wrfshtoi^'iKshwA"! d 
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3 
r 


This   is   also   a  two-center   integral.    The   result    is: 

\?        V  ,    HXu+v+a,-2v,      m     2u      £    -m-2v.      £,m        £2m 

L        h      x    (-D  2         2    Cu  Cv     V  wvl        wv2 

£2-m-2v2    ,    „,£   +£,-2v1-2v,-2u+v-s-r    _v+r+s 
x   C0  (kR)    1212  £ 


:{2(1-?2)    gJa2(£1-2v1-2u+v-r+3,£2-2v;d-s,C) 


^|  g3.    o    (£1-2v,-2u+v-r+2,£z-2v2-s,0}}    (4.13) 


where  the   g3-functions   are  defined  as: 

g.3    ,    (£>,£,, O    =12.     I         IL    (_l)u+v+ki+k2     (£1+|,2k1+u) 

X,l*<2  M=0    V=o     ic«©  ^s* 

x    (i2+*,2k2+v)    P^^'3)(£3+2k1+u,£lt+2k2+v,C) 


and  pJuJV,3\r,s,C)  are  defined  in  terms  of  the  sine  and 

X  iX,2 

cosine  integrals  as  below: 


»i"i;,s><^.-.»' 


=  {COSt.RCf^-TT   )(«U+V|0+«U+Vp2) 
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+  sin(KR£+i^Tr)<5u+V)1(6v0-6vl)}  J       _  cos  *„  ^  dt 


r+s 


+  {sin(KRC+-V^7T>(6u+v,0+6u+v,2) 


+{cos(<Rg+  \    2tt)(6  ■     -6      ) 
2    '   u+v,0   u+v,2 

+  Sln(KR5^_^>Vv,l,JM(W(i))    (2<RS-t)rts  dt 

(4)  IA0:  L  .'n(o«"I.(fi»t.(B0>,!1lJ0l"rtJ"i  ^ 


Here  we  have  a  two-center  integral.  We  want  to  ex- 

•2 

that  we  can  reduce  the  integral  into  sums  of  one-center  inte- 


press  j.  (<rn)Y   (rn)  in  terms  of  functions  centered  at  A  so 


grals.  Such  an  expansion  is  derived  in  Appendix  E: 


i.^'o^iV;'-1'1'4'!^'' 


-ju-r 


5  ' 
ICr'.nT'jJ^.m;*.1  ,m' )  (-1) 


m   -*■  in   "*■ 

Ji,(KRA)YJi,(RA)j£n«rA)Y£„(rA)      (4.14) 


for  r   in  the  intersphere  region 
With  the  above  equation,  we  have: 


-40    —   1V   ^    £   , 
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I^  =  (-l)S'24ir  V    T]   (-l)k  I(r\m;Ji2,m;S,',0) 


0  -».  r  m  .*.  m  ^. 

j£,(<RA)Y£,(RA)  j^  J£2«rA)Y£2(rA)Jr,«rA)Yr,(rA) 

COS0A   ?  • 

--T^  d3r  (4-15) 


Thus,  we  reduce  I._  to  an  infinite  sum  of  one-center 
integrals.  Since  j.,  decreases  rapidly  with  increasing  &', 
we  can  expect  that  the  sum  would  converge  rapidly. 

To  evaluate  the  integral  on  the  RHS  of  equation 
(4.15),  we  employ  the  same  technique  as  before  and  obtain 
the  following  expression: 


jii^i(KrA)Y^(?A)j^(KrA)Y^(?A)coseA/rA  d3r 
-N^<-l)£'+1CR/2>      WI     ■     E  I  I 

*i      X2  J|  ^o  Vjto  y=o  JsO 

x(-l)V'm  ^2m  Cm  (X+*2-2m-2v+2r+l    (kR)s(1_c2)S 

V  j .  V2  r        s 

£1+£z-2m+l-2v+2r-s    r  u  v  2        „       '  ^    ■ 

^  (    4C    BJiJlj£(8+3,C)    "-icKSll£lC«+2.e>H. 

(4.16) 
where  v  =   v   +  v2 
and  the   g**-f unctions   are  defined  as: 

V     t~      -c"»  u    (u,v)    (u,v,4) 

u1  o    v:e        mo 
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where  P  u:v'  ^(£,0  are  defined  in  terms  of  the  sine  and 
cosine  integrals  as  below: 


+  cosiC^-*,)*  6u+Vfl(«u0-«ui)>Jw   f(|)f^dt 


+  {cosi(£1-£2)TT  (6 n+6 rt) 


u+v,0  u+v,0' 


+  si 


n4(A1-H2)Tr  6      (6   -6   )}f      ^^^ 
u+v,l   ul  uO  J^-w)      t* 


dt 


■+{sini(£l+&2)n  (<SU+V)2-6U+V>0) 
+  cos£(Jt,+S,,)7T  6,,^  ,}        dt 


2J  U+V,1\L„  ,„Ky     t1 


(5)  TB0:   j   \1(Kr-)YT.^-)Jo.(<^)Y^(?n)cos0A/r^  d3r 


/ll  ~x        *   Byj£2    0'  £2   0'     A'  A 


rn  -*■ 
Expanding  j   (<r  )Y   (r  )  in  terms  of  functions  cen- 
tered at  at  B,  we  have, similar  to  IA0,  the  following  equation 

I*   =  (-l/247T  £   ]T   (_i)k  l(r<,m;£2,m;£',0) 

m  -*    r  m  ->  m  "*" 

x  Jr(<RB)Yr(RB)     nM(<rB)Y^(rB)j£„(KrB)Yr,(rB) 

2   3 

x  cos0A/r.  d  r 
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The.  integral  on  the  RHS  of  the  above  equation  can  be 
evaluated  as  below: 


I   VkVV?b>V(kVV(VcosV^  d3r 

hi     .  1  l  ■ 

=  ^  <  (-i)l.+1(R/2)  I  dC  I     E     E      Z 

*1+A2-2m-2v+3r  ^im        l2m  cm     ^^m^v-^r  a 

Vj  v2        r        s 

x(.l-5    )     e    l      2  .{(1+5    >8j|£j8(s+l,5)-  ^g|i&2(s,C)} 


where  v  =   v   +   v 
i         2 

and  the  g5-functions  are  defined  as: 

where  the  P     '   (£.£)  are  defined  in  terms  of  sine  and 
cosine  integrals  as: 


PiV'5)(£''5)  •'lsini(£*-^)ir  (6u+v,0+6u+v,2) 


+   cosi(lrl2)ir    «u+V)1(6u0-6ul)} 

f"^^      cos   t        _ 
x  I  2-  dt 

J  .■■I2icR?-t)    t 

♦  {cosiC^-l,)*    («u+V(p+6u+Vf2) 
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u+v,l      ul      uO 


1 


«(*+*<*>) 


sin   t  dt 


U(j-jo)/2<R^~t)  t 


+  {8ini(£1  +  fc2)ir    (6u+V)2-6u+v,0) 


+   cosi(£1+*2)ir    6u+vl}x  J  (2<R?_t)2p"  « 


(6>  '5o:  J„  v'"0^/*0"^"0^,'*0'00"6*'** " 

Expanding  j.(Kr0)Y£(r0) about   A,    we   have: 

i  T*  ■  V-  0  -»■  0       ->    v    l.  , £,  ,m 

^0  -     U    Js(KRA>Ys(RA)Jr-s(KRA)Yr-s(RA)Cs!r-l 


3r 


where     C^':2'"'     =    (-l)""1'^    (4tt)"  £  V 


*'>*3  rip'-Al     A&-U 

r+j.+jr-ar,   wa**& 


{    (-l)Pl+P2I(£",m;£i ,m;£' ,0) 


(m) 
IU^m^.m^j.O)    I    „         } 


I»  II  1*  H 


x   Y      (rA)cos0A/rA   d   r 
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Thus,  Iqq  is  reduced  to  a  sum  of  one-center  integrals 
which,  by  employing  the  Prolate  Spheroidal  Coordinates, can 
be  further  reduced  as  below: 

jjiJJti(KrA)j^(<rA)Y£3(?A)cos0A/rA  d  r 

=    (2£,  +  l)*    **    (R/8)  d?{    L        X       %3°   43"2V+1    (<R)r 

*<1-  sV  ^3"2v+l-r  {^  p^j^dt)^,^^^ 

*«<}-*<»)      tr+1 

,.  ,  r,x   F  J  P.   (§t)jn   (Jt) 

-(1/<R)  _±J ±-2 dt}> 

To  conclude  this  section,  we  note  that: 
(i)  We  have  reduced  the  three-dimensional  multi-center  inte- 
grals into  a  sum  of  one-dimensional  integrals  whose  inte- 
grands involve  the  sine  and  cosine  integrals.  The  sine  and 
cosine  integrals  can  be  evaluated  quickly  and  accurately  and 
can  be  called  directly  from  the  IBM  Scientific  Subroutines 
Package.  We  shall  discuss  such  evaluations  in  Appendix  G. 

(ii)  It  looks  as  though  the  expressions  for  I   ,  I   ,1   , 

s  AA   BB   AB 

I   ,  I   and  I   are  very  formidable,  yet  in  the  actual  cal- 
culations they  are  much  more  simple.  This  is  because  the  I- val- 
ues used  in  the  basis  of  the  MSXa  spin-orbitals  are  usually 
small,  and  in  most  cases,  do  not  exceed  3.  Thus,  for  example, 
if  we  use  s  and  p  atomic  orbitals  (I   =0,1)  for  the  a  state 
(m  =  0)  in  the  MSXa  calculation,  then  by  equation  (4.11), 
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there  is  only  one  term  in  the  summation  of  v,,  v2and  r  since 

a-m)/2  =0  .  . 
(iii)  By  partitioning  the  region  of  integration  into  differ- 
ent parts,  we  can  find  the  contributions  to  the  total  force 
by  the  electronic  charge  in  each  of  the  regions  and  for  each 
molecular  state: 


F      =   /   n  F  +  F 
total     4^   l     rN 


i 

=   "V  n  (F1  +  F1  +  F1  +  F1   )  +  F 
Z-   l^A    rB    rII    out;    rN 


where   i  represents  the  molecular  state  i, 

F.  is  the  force  exerted  by  the  electronic  charge 
inside  sphere  A  for  the  i   state.  Similarly  for  FR,  FTT 

and  Fout - 

F.T  is  the  force  exerted  by  the  other  nucleus. 

N 


4.5  Results  and  Discussion 
The  Hellmann-Feynman  force  calculations  have  been 
carried  out  for  two  different  systems:  the  hydrogen  molecule 
and  the  nitrogen  molecule.  The  hydrogen  molecule  is  investi- 
gated because  it  is  the  simplest  diatomic  molecule  and  accu- 
rate results  had  been  previously  obtained  (see  Kolos  and 
Wolniewicz  (1964))  so  that  we  can  make  some  comparison.  The 
nitrogen  molecule  is  investigated  because  its  electronic 
configuration  for  the  ground  state-  is  a  closed  shell  config- 
uration in  which  case  the  MSXa  method  should  work  better. 
Also,  it  involves  the  Tr-states  so  that  we  can  see  how  our 
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method  works  for  a  more  complicated  case. 
Hydrogen  molecule:  H2 

(A)  The  ground  state  for  the  hydrogen  molecule  H2 

1_  +  2 

is  l „   and  the  electronic  configuration  is  log.  Forces  for 

three  different  bond  distances  are  computed  in  order  to 
plot  a  force  curve  and  to  obtain  the  equilibrium  separation 
Re  and  the  force  constant  k.  First,  for  each  bond  distance, 
a  non-spin-polarized  MSXa  calculation  is  carried  out  by 
using  s  and  p  orbitals  in  the  atomic  regions  and  s  and  d 
orbitals  in  the  outer  sphere.  Then,  the  electrostatic  force 
is  calculated,  according  to  the  method  described  in  the  pre- 
vious section,  by  using  the  converged  wavef unctions.  The  re- 
sults are  listed  in  table  (4.1),  and  the  curve  is  plotted  in 
figure  (4.4). 

(B)  In  order  to  check  the  validity  of  the  Hellmann- 
Feynman  Theorem  of  the  MSXa  spin-orbitals,  the  non-muffin-tin 

total  energy  <E„(p)>  of  H0  is  calculated.  From  the  values 

•x,                             8<Exa(p)> 
of  <Exa(p)>,  we  obtain  the  values  of g-5- — —  at  the  three 

different  separations  by  least  squares  fit.  The  results  are 

listed  in  table  (4.2).  Two  such  curves  are  plotted,  one  for 

the  wavefunction  with  basis:  £=0 , 1  in  the  atomicsphere,  1=0,2 

in  the  outer  sphere;  another  for  the  wavefunction  with  basis: 

£=0  in  the  atomic  sphere,  £=0,2  in  the  outer  sphere.  Also, 

to  compare  with  the  experimental  result,  a  curve  with   -  — 

against  R  where  E  is  the  total  energy  calculated  by  Kolos 

and  Wolniewicz  (see  kolos  and  Wolniewicz  (1964))  is  plotted. 

Such  an  H2  calculation  is  considered  to  be  more  accurate 
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TABLE  4.1 

Force  Calculations  For  H 

2  <Yl 

2 

■2F?g. 
A 

2F°S 

2F0g+ 
out 

2Fjf 

FN 

F 
total 

R=l 

4 

0.264 

0.271 

-0.026 

0.442 

-1.020 

-0.068 

R=l 

44 

0.262 

0.266 

-0.U26 

0.417 

-U.965 

-U.046 

R=l 

5 

0.257 

0.258 

-0.025 

0.403 

-0.889 

0.004 

*  Basis:  £=0,1  in  the  atomic  spheres 

£=0,2  in  the  outer  sphere. 
a=0. 77627 

R  in  atomic  unit;  F  in  Rydberg/a.u. 
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TABLE  4.2 


3<Exa(p)> 


total 

3R 

l2 

(a) 

F 
total 

a<ExaU) 
8R 

(b)                  %(c) 
>      _9<EXa(p)> 

dR 

(d) 
3<E> 
8R 

R=1.4 

-0.068 

-0.024 

-0.035 

-0.004 

R=1.44 

-0.046 

0.004 

-0.01^ 

0.032 

R=1.5 

0.004 

0.046 

U.024 

0.068 

Re 

1.4y5 

1.43 

1.46 

1.401 

ke 

0.72 

U.71 

0.60 

0.73 

(a)  Hellmann-Feynman  Force  with  basis:  1=0,1   in  the  atomic 
region;  1=0,2   in  the  outer  sphere, 

(b)  Basis:  &=0,1  in  the  atomic  region,  £=0,2  in  the  outer 
sphere. 

(,c)  Basis:  £=0  in  the  atomic  region;  J6=0,2  in  the  outer 

sphere, 
(d)  Calculated  with  the  Kolos-Wolnlewicz  wavefunction 
Force  in  Ry/a.u.;  K  in  atomic  units;ke  in  Ky/a.u. 


Fig.  4.4   Force  Curves  for  H2 

1.  Basis:  X,=0,1  in  the  H  sphere 

«,=u,2  in  the  outer  sphere 

2.  Basis:  £=0  in  the  H  sphere 

1=0, -A   in  the  outer  sphere. 


fc>7 


Force(,Ryd/a.u.  ) 


0.U6     [_        X    Experimental 
n     aExa(p)       * 
9R 

aEXct(£)     2- 

9R 

0-04   L     a  Pjm(?) 


U.02 


-U.02 


-0.04 


-0.U6 


-U.08 


-0.10 


1.2 


1.3 


1.4 


1.5    separation  (-a /u.  ) 
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TABLE   4.3 


£-D'ependence  of  Hellmann-Feynman   Forces    for  H2   at   R=1.4   a.u, 


„     aS"  afr  acr  aS 

2F/      2FBS      ,FQgt     2Fir      FN      Ftotal 


ffl 

0.265 

0.271 

-0.026 

0.442 

-1.U20 

-0.068 

*2  • 

0.268 

0.274 

-0.026 

0.442 

-1.U20 

-0.062 

#1:   £=0,1  in  the  atomic  region;  1=0,2    in  the  outer  sphere 
#2:   £=0,1,2  in  the  atomic  region;  £=0,2  in  the  outer  sphere. 
Force  in  Ry/a/u. 
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than  that  experimentally  observed.  The  equilibrium  separa- 
tions Re  and  the  force  constants  ke  for  each  of  the  curves 
are  obtained.  The  results  are  listed  in  table  (4.2)  and  the 
curves  are  plotted  in  figure  (4.4). 

(C)  To  check  the  sensitivity  of  the  force  with  re- 
spect to  the  basis  functions,  two  different  calculations,  one 
with  s,  p  orbitals  in  the  atomic  region  while  the  other  with 
with  s,  p,  d  orbitals  (both  with  s,  d  functions  in  the  outer 
sphere),  are  carried  at  R=1.4  a.u.   The  results  are  in  table 
(4.3). 
Nitrogen  Molecule:  N„ 

Bader .Henneker  and  Cade- have  done  some  force  calcula- 
tions by  using  Hartree-Fock  wavefunctions  as  the  basis  (see 
Bader,  Henneker  and  Cade  (1967)).  In  order  to  compare  these 
two  methods,  a  force  calculation  for  nitrogen  molecule  is 
also  done  by  using  the  MSXa  wavefunctions. 

Both  the  MSXa  and  the  Hartree-Fock  calculations  are 

carried  out  at  experimental  equilibrium  separation  R  =2.07 

1  + 
a.u.   The  ground  state  for  N2  is   Z   and  the  electronic  con- 

2   2   2   2   4   2 
figuration  is  la  la  2a  2a  Itt  3c  .  The  spin-  polarized  a  value 

g   U   g   U   U   g 

for  nitrogen,  q=0. 74522  is  used,  and  the  I    values  used  for 

the  different  molecular  orbitals  are  listed  as  in  table  (4.4). 

Here  we  treat  the  la„  and  la  as  core  states.  The  results  of 

S  u 

the  force  calculations  are  as  in  table(4.5). 
Discussion  : 


When  we  compare  ^(p),  Fx  (p)  and  the  experimental 
values  for  the  force,  we  find  that  we  have  achieved  a  big 
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TABLE  4.4 


£-Vaiues  in  the  MSXct  Calculation  for  N„  (  £  ) 


2a 
g 

2a 
u 

lTT 
U 

3a 
g 

atomic 

"0,1,2 

0,1,2 

1,2 

0,1,2 

outer 

0,2 

1,3 

1,3 

0,2 

*  la   and  la   are  both  treated  as  core  states. 
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TABLE  4.5. 


Hellmann-Feynman  Forces  for  N2  at  R=2  a.u. 
Calculated  with  MSXa  and  Hartree-Fock  Wavefunctions 


2   2 

la  la 
g  u 

2 

2a 
g 

2 

20 
u 

4 

ITT 
U 

2 

3a 
g 

* 
fA 

0. 

1.513 

-0.501 

1.278 

0.462 

fB 

7. 

2.227 

0.9b2 

2.317 

l.b85 

f   + 
out 

0. 

U.0U3 

-0.844 

0.160 

-0.788 

fII 

0. 

3.3±6 

-0.0U3 

3.274 

0.452 

xa 

7. 

7.059 

-0.395 

7.029 

1.811 

fHar- 

-Fock 

7. 

858 

9.387 

-1.S21 

8 .  512 

0.525 

Ftotal(pxa 

)  = 

-1 

996 

Ftotai (Hartree-Fock)  =  -0. 

*  Total  electronic  force 
Force  in  Ry/a.u. 
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improvement  by  using  the  non-muffin-tin  charge  density  to 

calculate  the  force.  The  force  constant  k  and  the  equili- 

e 

brium  separation  for  H  calculated  with  this  method  involve 
errors  of  1.4%  and  7%  respectively,  while  the  FXct(P)  indi- 
cates no  binding  at  all. 

'v        9Exct(p) 

However,  when  we  compare  F._  (p)  and  - —  ,  we 

Aa  9R 

find  that  there  is  some  difference  between  them.  This  means 
that  the  Hellmann-Feynman  Theorem  is  not  quite  correct  with 
the  approximate  MSXa  wavef unctions.  Some  reasons  may  be : 

(1)  Since  the  wavef unctions  we  used  are  hot  exact 
solutions  to  the  Xa  Schrodinger  equations.  It  also  assumes 
the  muffin-tin  potential  approximation  and  the  potential  is 
generated  only  from  a  muffin-tin  averaged  charge  density. 
Thus,  it  is  only  an  approximate  Xa  wavef unction  and  we  cannot 
claim  that  it  should  satisfy  the  Hellman-Feynman  Theorem  very 
well.  Actually,  due  to  the  polarization  effect,  the  potential 
on  both  sides  of  the  atomic  spheres  should  be  different.  The 
potential  in  the  bonding  region  facing  the  other  nucleus 
should  be  less  shallow  than  the  potential  in  the  outer  re- 
ion.  Thus,  there  should  be  more  charge  in  the  bonding  region 
than  what  is  expected  by  using  the  muffin-tin  averaged  poten- 
tial. This  means  that  the  actual  force  should  be  more  attrac- 
tive than  the  force  we  got  by  assuming  the  muffin-tin  poten- 
tial. This  polarization  effect  should  be  especially  severe 
for  the  core  state,  for  which  the  electronic  charge  accumu- 
lates near  the  nucleus. 
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(2)  Besides  assuming  the  muffin-tin  averaged  poten- 
tial, we  also  use  only  a  finite  set  of  basis  functions  with 
small  £  values.  This  may  not  be  important  in  calculating 
the  total  energy  but  it  may  be  important  in  calculating  the 
force,  because  the  force  is  much  more  sensitive  to  the  wave- 
function  than  the  total  energy  is.  We  further  note  that: 

(i)  In  the  force  calculation  for  the  hydrogen  mole- 
cule  (  table  (4.2)  and  figure  (4.4) ) ,  the  values  of  F   (p) 
by  using  the  basis  (£tt=0,1;  £   t=0,2)are  closer  to  the  values 
of ■  3^xal£2  calculated  with  the  basis  (£R=0;  *o'ut=0,2)  than 
to  that  obtained  with  the  basis  (£H=0,1;  £out.=0,2).  Further- 
more, we  have  the  following  relation: 


'w^  (V°'1;  Vt=°.2>  <  %—  <V°;  'out'0'2' 

<  9E^ip)  (£„=0,1;  £   +=0,2)  <   F         ■■  . 
3R     H       out  experimental 


Or 

This  can  be  expected  since  if  the  error  in  F   (p)  (£H=0,1; 
£ou-t=0,2)  is  of  first  order  of  the  errror  in  the  wavefunction , 

then  the  error  in  E   (p)  (£R=0,1;  £Out=0-2)  is  of  second 

>\, 
order  and  that  the  error  in  E   (p)  (£=0;  £   =0,2)  should 

xorK/   H     out 

be  larger  than  that  of  E   (P-)(£H=0,1;  ^out=0>2)  but  smaller 

than  second  order. 

(ii)  In  table  (4.3),  we  see  that  when  we  increase 

the  i -values,  the  net  force  becomes  less  repulsive  and  agrees 

better  with   xqiPJ  #  The  improvement  is  very  slight.  However 

3R 
if  we  use  only  s  orbitals  in  the  atomic  spheres, then  F^=0 
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and  the  net  force  becomes  more  repulsive.  Thus,  the  p  orbital 
is  very  important  in  the  force  calculations,  although  the  s 
orbital  is  almost  sufficient  for  the  total  energy  calculation. 

In  the  force  calculation  for  the  nitrogen  molecule, 
we  see  that  the  MSXa  forces  have  the  same  characteristic  as 
the  Hartree-Fock  forces.  Again,  the  MSXa  forces  are  too  re- 
pulsive. However,  the  relative  values  of  forces  contributed 
by  the  different  molecular  states  for  both  cases  are  compar- 
able, and  the  au  state  in  both  calculations  contributes  re- 
pulsive force.  Of  course,  we  should  not  expect  our  result 
could  be  the  same  as  that  of  the  Hartree-Fock  calculation. 
What  we  should  expect  is  just  a  general  characteristic  of 
these  two  kinds  of  forces. 

In  conclusion,  the  use  of  MSXa  wavef unctions  instead 
of  the  muffin-tin  averaged  charge  density  to  calculate  the 
force  makes  a  big  improvement.  It  gives  a  reasonable  equili- 
brium separation  and  a  good  force  constant.  However,  due  to 
the  muffin-tin  potential  and  the  truncation  of  ^-values  in 
the  basis  set,  the  existing  calculation  of  MSXa  wavefunctions 
still  has  to  be  modified  to  satisfy  the  Hellmann-Feynman 
Theorem. 


CHAPTER  V 

D I POLE  MOMENT  CALCULATIONS 

5.1  Introduction 
In  Chapter  III, we  have  discussed  how  to  compute  the 
dipole  moment  of  a  diatomic  molecule  by  assuming  the  muffin- 
tin  charge  density.  In  this  chapter,  we  shall  discuss  the 
same  calculation  but  using  the  non-muffin-tin  charge  density, 
i.e.,  the  charge  density  formed  from  the  MSXa  wavefunctions . 
We  shall  first  discuss  the  computational  procedure  and  then 
give  the  results  for  some  diatomic  molecules.  The  results 
show  that  the  use  of  non-muffin-tin  charge  density  indeed 
provides  big  improvements  over  the  use  of  muffin-tin  charge 
density. 

5.2  Method  of  Dipole  Moment  Calculations 

To  calculate  the  dipole  moment  of  a  diatomic  mole- 

t  '     -*■  3 

cule,  we  need  to  evaluate  the  integral    p(r)z  d  r   contrib- 
uted by  the  electronic  cloud.  Again,  in  the  MSXa  calcula- 
tions, we  have: 

p(r  )  =   2_  ni  ui(r)  u±(r) 
I 
As  usual,  we  partition  the  coordinate  space  into  three  re- 
gions: the  atomic  region,  the  intersphere  region,  and  the 
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out 


Fig.  5.1   Schematic  Representation  of  A  Heteronuclear 
Diatomic  Molecule  in  the  MSXa  Calculation. 
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outer  sphere  region  (figure  5.1).  As  we  have  mentioned  in 
Chapter  III,  the  value  of  the  dipole  moment  is  independent 
of  the  choice  of  the  origin  for  a  neutral  molecule,  we  can 
choose  the  center  of  the  outer  sphere  as  the  origin.  Then 
we  have : 

.jp(r)z0  d  r  =  r  n±    (u*  +  y£  +  yout  +  yII}        (5'1) 

i     f   %*  ->  %  -*-    .  ; 
where        yA  =  J   ui(r)ui(r)z0  d 

yB  =  J,  Vr)ui(r)zO  d 

yout  =  JoutUi(r)Ui(r)Z0  d  r 
WII  =  JrI  ui(r)ui(r)zQ  d* 


3 
r 


,3 

r 


.3 

r 


and  the  dipole  moment  u  is:      u  =  ■  u.T  -  u 

N    e 


where  mn  =  ZgRgQ  -  Z.R.q   is  the  nuclear  contribution, 

f   ■*■  3 

y„  =    P(r  )Zq  d  r   is  the  electronic  contribution, 

Z^,  Zg  are  respectively  the  atomic  numbers  of  the  two 
nuclei  A  and  B;  u  is  in  atomic  units. 

Integral  Over  Sphere  A:  u 

Inside  sphere  A,  we  have: 
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u.(r)  =  /_.  C   R  (r.)Y  (r„) 
1       X   Am   £   AA  JT  Ay 


v,..,->  r  Y°(?  ,  -  „ 
0    A    AO    ^"/-/   rA  luAy     AO 


and     zn   =  z.  -'R.A  =  (4tt/3)2  rAYj(rA)  -  R 


Thus ,  we  have 


,3 

r 


yA  =   JA  Ui(r  }  Ui(r)  Z0  d 

-   I  ?  Cifm  <m  j  R£1(rA)Y;i(?A)R£2(rA)Y;2(?A) 


x  {  (4T/3)*rAY;(?  )  -  RAQ}  d3r 


"\'\    C^m  C^m  {  (4*/3)*  IC^.m^.mil.O) 


x  (  \z    (rA)R„_(rA)r?  dr, 

Jo 


A'"A2V.  AyiA  "XA 


1*2  J  R£, 


"  RA0  h,l7         Ri    <rA>R£  <rA>rA  drA  } 


Integral  Over  Sphere  B:  yR 

Inside  sphere  B,  we  have 


V^  '     £  C«  Ri(rB)  .Tj(fB) 


and     z0  =  zB  +  RB0  =  (47T/3)*  rB  yJ(?b)  +  RB0 


Thus ,  we  have : 
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i  f      %*  -»■     ^     -*■ 

U_   =  u.(r)    u.(r)    z      d 

B         JA        1  1  0 


3 
r 


t->    r-        iB         iB  a 

=      >     >      C  C  (4tt/3)2    l(£1>m;£,,m;l,0) 


j     j        S,  j  m      2, ,  m 


rRn  Q  r*> 


3   .  _  r»      ,     v„     ,     .2 


X 


rV^V^B    drB   +   RB0fi4l£2    f  \ 


(rB)R£2(rB)rBdrB> 


Integral  Over  the  Outer  Sphere:  uout 

Inside  the  outer  sphere,  we  have 


^   -»■     V   10  m  -»■ 

V*>  "  £  c*  .  YV  Y?(ro> 

Thus      u      ^   =  u.(r)u.(r)z_d 

out  Jout     11  0 


,3 
r 


x    (47r/3)*r0Y°(?0)    d3r 
■     |    ^   Ci?mCi°n   (47T/3)      K*».n,;£afm;l,0) 

XT     V^W'o  dro 

Integral  Over  the  Intersphere  Region:  ujj 

In  the  intersphere  region,  we  can  express  u.(r)  as 
linear  combinations  of  products  of  Neumann  functions  and 
spherical  harmonics  as  discussed  in  the  force  calculations. 
In  the  same  way  as  before,  we  obtain: 
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,3 
r 


yil   =  Jxi  Ui(r)  Ui(r)  Z0  d 

+  V  aJ°  A„°   IrtA(*,,a2,nO  ■+  2  X  AlA  AlB  IAO(a,,a,,m) 

+2  £  A^A  A*°   IA_(£,,A  ,m)  +  .2  2  A*8  A*°  IDrt( *- ,  .  *,  .m) 
t   i      £  n  £  ni  AO   i'  2'       r-'i  8,  -m  £  m  BO   i   z' 

(5. J 
where  I^Ci,.^  .m)  -  .]■    n£  ,  (^Y^  (r^n^  (-^A^ 


(5.2) 
,3 


V'.'V"  "  i„ni1<KrB>^l(fB)nt1<"B)TIt<fB>V3' 


IABU,,t,,m>  -  [^  nJl.(KrA)Y7i<ifA)ni2(KrB)^(fB)z0d 


3 
r 


3 

r 


IB0(^1^2.m)  =  J  n£  (KrD)Yf  (r„)jn  (<rrt)Y™  (r„)znd3r 


•b-1^V^2(«Vy?2(Vzo^ 


To  evaluate  the  above  integrals,  again  we  employ  the 
Prolate  Spheroidal  Coordinates  as  we  did  in  the  force  calcula- 
tions. However,  since  in  this  case  we  have  a  heteronuclear 
diatomic  case,  so  the  limits  of  integration  when  transformed  to 
Prolate   Spheroidal   Coordinates  would  become  more  complex. 
We  shall  discuss  each  of  the  integrals  in  the  following: 


(1)    I      :       f       n      (<rfl)Y1J1   (rA)n,    (*rA)Y™   (rA)zft  d 
AA         JIX       x,]         A      *!      A      x,2         A      x,2A      0 
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3 

r 


This  is  a  two-center  integral .  We  shall  express  the 
integrand  in  terms  of  Prolate  Spheroidal  Coordinates  with  A 
and  0  as  the  two  foci  (f igure(5. 1) ) : 


?  "  (rA  +  r0)/RA0 


n  =  (rA  -  r0)/RA0 


JA        *0 


The  reciprocal  relations  are: 


rA  =  iRA0  (?  +  Tl)  ■■r0-*?A0(*-n) 


cosGA  =   (l+Sn)/(S+n)  and  cosO     =  -  (l-£n)/(S-n) 

sin0A  =  «g  zlHU  )>  sin9o  =     US  -Dd-n  »* 

(5.3) 
and  the  volume  element  is: 


d3r  =  (i  RAQ)3  (£2  -  n2)  d£  dn  d<}> 


Since  the  two  atomic  spheres  are  in  general  of  dif- 
ferent size,  the  geometry  of  the  partition  is  a  function  of 
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two  parameters,  the  bond  distance  R  and  the  ratio  of  the 
radii  of  any  two  spheres,  in  this  case,  we  choose   6  =  R3/R. 
Then  we  can  express  the  different  quantities  in  terms  of 
these  two  parameters : 

RA  =   (l-3)/R:   RB  =  8R;   RQut  =  R;   RAQ  =  BR; 

RB0  =  (l-6)/R  <5.4) 

With  the  coordinate  system  as  expressed  in  equation  (5.3) 
and  also  by  the  relations  in  (5.4),  we  have  the  following 
restrictions  for  the  intersphere  region: 


(1)   rA  >  RA 


implies   £RA0(£+n)  >,     RA      or  £+n  >  2(1-B)/I 


<2>   r0  *  Rout 


implies   §RAQ(S+n)  S  Rout     or   C-n  *  2/6 


(3)   rB  >,RB 


2     2     2 

and  r_,  =  r,  +  R  -  2r.Rcos9. 
B  •  A  •         A     A 


implies    n  -  1  (1-  i6)Cn  +  S  -  %   (6  +3-1)  >,     0 


(4)   £  >,   1 


(5)   -1  $  n  *  1 
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~»£ 


1  - 


-l  - 


(ii)  6<l/2 


L-M 


Fig.  5.2   Boundary  of  the  Intersphere  Kegion  in  Prolate 

Spheroidal  Coordinates  with  A  and  0  as  the  Foci. 
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Combining  (1),  (2),  (3),  (4)  and. (5),  we  have  the 
transformed  region  as  defined  by  the  following  limits  (figure 
5.2): 


5  *  5  = 

o 


L?-  3 


if   6  >, 


if   6  < 


5  *  T  +  1 


f  (C)  $  n  *  f  (■£) 

1  2 


where    f  (?)  =   r{(1-46.)5-  {(l-B)C  +  (B  +6-1)}  } 

2  P 


r  2/6  -  2  -  £     if   5  <  2/3-1 


and     ■  f  (?)  = 

i 


L  5  -  2/| 


if   £  >  2/6-1 


Expanding  the  spherical  harmonics  and  the  Neumann 
functions  in  terms  of  £,  n ,  <t>  as  we  did  in  the  force  calcu- 
lation in  Chapter  IV,  we  have  the  following  result  after 
some  straightforward  algebra: 


V^=o 


Wiw-aV**'- 


5"   .  ,    r   J^m   A2m  _m   i1+&2-2m-2v+2r      . 
£_j   i(-l)   u    w .    C  C  (  KA0; 


s-2 


V,     V. 


r  s 
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*  1  ^2 

2 

-  ^— —  h^£.(s,0  +  C/(kRaq)2  hJ^Cs-l.C)}} 


RA0 


where  v  =  vx+  v2 


and  the  hx-functions  are  defined  as: 


U+V   (u,v) 
An 


QjUIV'1)(s+2n+u+v)0 


where  the  A-coef f icients  are  also  defined  in  Chapter  IV  and  the 
Q-f unctions  are  defined  in  terms  of  sine  and  cosine  integrals, 
similarly  to  the  corresponding  terms  in  the  force  calcula- 
tions except  for  the  integration  limits: 


Qi"il,1)(1»5)  =  (  «>'«*<  W  <au+v,0-6u+vf2) 


KRCVt^")) 


\  cost 

+  sinf(£  +  £  )tt  5U+V  !>  I  -rr- 

1    2      U+V,X   JtR^+f(p)       tx 

+    {s.ini(A1+£2)ir    (5u+V)2-6u+v0) 

+   cosi(£.+£,)TT   6   ^      .}  sint 

1       z  u+v    II  


dt 


+    {sin|(£1-£2)Tr    5„_r   1(6„n-6ir1) 


u+v,lv"vO~  vl 


dt 


tdCJ+MO 


/">>«w/    d 

+  cosi(£1-£2)Tr    (6  6  )}  - 

u+v,u     u+v,z  W^f.cy)    t 


dt 
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(2)  IAn:  jn  n.  (<rfl)Y:  (?A)j0  (Krn)Y™  (?n)zAd^ 


m  ,-»•  N  .   ,    ,„ni  / ->  x   ,3 
AO-  Jn  ''£/^A-£/rA)jJl2(Kr0)Y*2(r0)20C 

This  is  also  a  two-center  integrals,  we  can  use  the 


same  coordinate  system  as  in  the  evaluation  of  IAA.  The  re- 
sult is: 

m   m      I  4-3  ft*'     ^  >  ^"^ 

IAA  =  N0   N0   (-1)  »  (RA0/2)  (<RA0)   J    d£   2.  L       Z. 


LAA  "  %  iN£2  .^-■L>    ^AO'^  ^RnA0 


}o  U«©   v  so    y>  —  o 
I.-M-XV,    litf,,-"*)]   V""***. 

f*      2-    {(-1)     03  »   .«  2  (-1)  2    2 

,v=o      Vv»«>       f^o  Vl      V2 

m  2u   £1-m-2v,   £,-m-2v2-      £ . -2v  -2u+v-r+4.-2v  -s 
x  C   C    C  :      *  C  2      2(<R.n)  -1    *  22 

u    v   r         s       v   ACr 


x  (^^-^i-^a^v.-s   cv+r+s  {(1+^) 

x  hJjJi2(Ji1-2v1-2u+v-r)£2-2v2-s,C)  -  5/-(kRA0) 
x  hl    I    (^1-2v1-2u+v-r-l,^2-2v2-s,^)}} 


where   h2  -  U,  ,£„,£) 


u+k,+k„ 
{  (-1)    *  -  2(£1+i,2k1+u)(Ji2+§)2k2+v) 


x  Q(tu'sv'2\l3+2k1+\i,i^+2k2+v,E>)    } 


and  the  Q-functions  are  defined  by: 
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Qp(U;V'2)(p,q,0  =  (sinCcR.^-iU^)*)  (6_  ~6  o) 


+  cos(KRA0C-i(^1-£2)Tr)  6u+vl} 

x  rw+^      cost     dt 


+  tcos(KRA05-*(£1-£2)ir)(«u+Vi2-«u+Vi0) 

+  sin(KRA05-i(a1-a2)Tr).  5U+V)1> 
[  sint       dt 


+  {sin(KRA0?+i(£1+£2)1r)  (6^  0+«u+v>2> 


+  cos(KRA0C+i(£,+£„)TT)  6 .,($__«-$„..)} 


1   2'    u+v,l  uO  ul 


rKfcof^^O 


X   )  ~S n  dt 

H,nH<rV  tp(2<RA05-t)q 

<3>  ZBB:  L^1^rB<1^B)^2(KrB)Y?2(?B)ZOd 


,   3r 

This  is  also  a  two-center  integral.  We  treat  this 
integral  by  considering  the  Prolate  Spheroidal  Coordinates 
with  0  and  B  as  the  two  foci: 

5  =  <r0  +  rB>/RBO  '   n  =  (r0  "  V/^O.  '  *  =  *0  =  *B  . 


With  this  coordinate  system,  we  have  the  following  restric- 


88 
tion  for  the  intersphere  region: 


<X>   rB  *   RB 


implies  (FL.n/2)(£  -  n)  *  Rn  or  £  -   n  £  28/(1  -  0) 


'BO'  /v^   .  "  '   B 


(2)   rn  $   R   . 
v  '    0  ^   out 


implies  (RRn/2)(C  +  n)'S  RA,1+.   or  .S  +  n  $  2/(1  -  B) 


W'MS  r  "'  *  "out 


(3   rA  >  RA 


2     2     2 

and  r^  =  rg  +  R  +  2rgR  cosQg 


24      Y        2   4  2 

implies    n  +  -  (1  -  i)  Sn  +  £  +  -  (1  -  y  -  Y  )  >  0 


where    y  =  1-3  =  R4/R 


(4)   C  £  1 


(5)   -1  $  n  S  1 


Combining  all  the  above  restrictions,  we  have  the 
transformed  region  as  defined  by  the  following  (figure  5.3) 
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(i)  Y<l/2   (0>l/2) 


Ui)  Y>  1/2  (g<l/2) 


-*t 


Fig.  5.3  Boundary  of  the  Intersphere  Region  in 
Prolate  Spheroidal  Coordinates  with  0 
arid  B  as  the  Foci. 
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5  >,  V   = 
o 


2/Y  -  3      if  y  <  i  (6  >  -i) 


if  Y  >  |  (8  <  |) 


C  <  2/y  +  1 


g(0  *  n  $  g,(£) 

1  2 


2_*. 


where   gj(C)  =  -  {-  (1-y/2)C  +  ((1-yK   -  (1-Y-Y  ))2> 


g  (C)  = 

2 


C  -  2/y  +  2       if   £  <  2/Y-  1 


2/Y  -  € 


if    C  >  2/Y-  1 


Expanding  the  spherical  harmonics  and  the  Neumann 
functions  in  £,,    r\,    <j> ,  we  have  the  following  expression  for 
IBB: 


XBB  -  "?,  <2    (-D*'+*'+i'  (HB0/2)4  f;  d5  (  J  'f "' 
y  V        (_1)£1+il2-2m+3r-2v   £xm   £2m 


Vl       V2 


x  cm  cAi+i2-2m+2r-2v        s-2   £ ,+A2-2m+2r-2v-s 
r  s  ^<nB0;     ^ 


x  (1-C2)S  {2CU+S2)  hf  o  (a+l.C)  -  ^V1  h*  £  (s^> 

1  2  R^«    i  2 


BO 


+  C/C^q)2  h»i£a(s-i,o» 
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where  v  =  \>1   +  v2 

and  the  h3 -functions  are  defined  as: 


u+v  >(u,v)^(u,v,3) 

uao  vso     no  0  "        ~  1   2 


,        ■  V*  ^T   v       u+v  Cu,v;  (u,v,3) 
*Ll ,(S'?)  "III    ("^    An    °<lL      ;(s+2n+u+v,0 


where  the  A-coef f icients  are  defined  as  before  and  the  Q- 
f unctions  are  defined  as  below: 

<C3)(£'°  =  {  C°S  *<  W  (au+v.0-6u+vi2) 

r«WJ-fc<)»  cost 


int  dt 


+  {sin*<*l+£a)ir  (6U+V)2-SU+V>0> 

■3-3 

+   <cosi<V*2)ir    («U+Vf0+6U+Vf2) 


sin 


dt 

7T 


- 

•<W3-3.oo  z 

(4)  W  (;i1(KrB)Y;i(rB)jl!(Kr0)Y;a(r0)z0d3r 


We  evaluate  this  kind  of  integrals  by  using  the  same 
coordinate  system  as  defined  in  the  treatment  of  IBB,  i.e., 
with  0  and  B  as  the  two  foci.  The  result  after  some  straight- 
forward algebra  is: 
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IB0  -  N?1   N?  C-l)£l(4R  n rW  )-3J       «  J    £      t 

BO           *i      *2  BO             BO         Jj;            ^0   ^      £- 

I   E  z    t(-i)u+vt*,-2v»-;!B  «;;n  cm  c2 

V°        r7%  fr*                                              vi     .    v2        u        v 


x   cix-m-2v1    2  )r+s-2u+v+2m  2.m+r+s   £v+£ !-m-2v x-r 


x   ?£2-m-2v2    s   {(1_c2}    hj^^g^u+v+m.r+m.C)    +    (5/^0 ) 


x   h„        (s-2u+v+m-l,r+m,C)}} 


where  the  h"1 -functions  are  defined  as: 

hj    .    (r,s,C)   -    III       L    (-D        *      2(Ai+i.2k1+u) 

x    (*2+J,2k2+v)   QjU;V'4)(r+2k2+v,s+2k+u,C) 

J6  iX,2 

and  the  Q-f unctions  are  defined  as: 


<X'4)(P'q'C)  =  {8in(,cRB0€-*(li+£»),r)(6u+v,0-au+vi2) 


+  cosCKRgQC-iC^i+^z)^)    5u+v  1) 
*  COSt/{(2KR      £-t)Ptq}    dt 

vg)-3.(v)  B0 

■+   {cos(<RB0C   -m+L^-n)    C6u+v>2-«u+Vf0) 


y3 


x  ,  .     sint/{(2<R      £-t)Ptq}    dt 


B0  *      2  u+v,0     u+v,2y 


x  ,  l/{(2KRnr.C-t)Ptq}    dt 


(5)  W  ju  h^o^^o^h^oK^o^o  d 


3 
r 


Since       ^   (?0)Yj   (?0)   =    £    1(1,  ,m;*2  ,m;*3  ,m3  >y£(?   > 
therefore      IQ0  =    £     K^  ,m;  £2  ,m;£3  ,m3 ) 

X    JIXJi/Kr0>J£,^r0)YTS/'0)a0   d3r 

Since  the   integrand  depends  on   <(>   only  through  Y^3(rQ) 
therefore,    after   integrating  over  <J> ,    the  only  nonvanishing 
term  is  the  one  with  m3=0.    Thus,    we  have: 


:00  =     Z      I(Ai,m;&2,m;£9fO)    f  j.    (icr    )j      (icr    )Y°   (r    )z.   d 
Vl'.-J.l  Jri      *  2  3 


3r 


Also,  since  Y°3<?0)z0  =  Y£3<?0)r0coseo 


-  (4it/3)V?3(v  Y°(v 
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=  (4ir/3)*r   T,      I(*3.0;1,0;^,0)Y°  (r  ) 
0  ^-W,-'l  *-   ° 

therefore   I   =  (4tt/3)2  £   .  £  I(  J.  x  ,m;  l2  ,m;  i3 , 0) 


x  I(£ii0;li0;*,f0)JjJli(Kr0)JJla(Kr0)Y^(?0) 


x  r  dJr 

Employing  the  Prolate  Spheroidal  Coordinates  as  in 
the  evaluation  of  IB0,  and  expanding  the  integrand  in  terms 
of  these  coordinates,  we  have: 


f   Js  (<rn)j0  (<rn)rnY^  (fft)  dur 


r0  ,£  x  ^3. 

'£1x"*0/u£2x'%*0/*0*Altvr0' 


0      h  4-3  c*^   £.0   'VJ/        r 

N^  (2lT)$(RB()/2)  «RB0)    £   ^     X!   (,CRB0) 

ys.o  f=  0 

j',   dC  {?^"2v"r(l-?2)r  {2C  J  "    J   (4t)JB  (4t)/tr~2  dt 

*•  vjin*)  *.        2 

(i/<RRn)l       ja   (4t)d0  (it)/tr-3  dt}} 


BO 


(6)  IAB:  )  n,i«rA)Ymi(?A)n,2«rB)Ym2(?B)Zo  d3r 

'XI 

This  is  a  three-center  integral.  But  we  can  write  zQ 
as  za-Ra0  and  reduce  it  to  a  two-center  integral.  We  shall 
evaluate  this  integral  by  using  the  Prolate  Spheroidal  Coor- 
dinates, in  this  case,  with  A  and  B  as  the  two  foci: 
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C  =  (rA+rB)/R  ;      n  =  (rA"rB)/R  ;       4>  =  <t>A  =  *B 

With  this  coordinate  system,  we  have  the  following 
restrictions  for  the  transformed  intersphere  region: 


(i)  rA  *  RA  °r     5  +  n  »  2(1  -  B) 


(ii)  rB  >,   RB  or     E,   -  n  >  23 


(iii)  rQ  $  R 


Since   rg  =  rJ  +  r2q  _  ^  R^   cosqa  ^   R! 


Therefore   ?2  -  (43-2Kn  +  n2  -  (4+43-432)  ^   0 


(iv)   ?  *  1 
(v)  -l  <  n  <  i 

Combining  the  above  five  restrictions,  we  have  the 
transformed  region  as  defined  by  the  following  limits  (fig- 
ure (5.4)): 


1  «  5  ^max  :     y/5)  *  n  *  7,(5) 
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(i)  8>l/2 


|(26-l)(i-"-)2l<l 


n=y,U) 


(iii)  B<i/2 


and  1(2B-D^1"BJ  >>\>1 


Fig.  5.4  Boundary  of  the  Intersphere  Region  in  Prolate 
Spheroidal  Coordinates  with  A  and  B  as  the  Foci 
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where 


'max 


,2.,-i 


2.  ,4 


r  {(1+B-B  )/(B-B  )f2     if  |(2B-l)l(l-B-B  )/(B-B  )}"<  1 


1  +  21 


I  3  -  2f 


if  B>|  and  (2B-1){(1-B-B2)/(B-B2)}*>  1 


otherwise 


y,(C)  - 


2(1-B)-C 


if  SS3-2B 


(2B-l)C-2{l-(B-B2)(^2-l)}i 


otherwise 


y  (5)  = 

2 


C-2B 


if  SS  2B+1 


L  (2B-D^+2{1-(B-B2)(C2-1)}2 


otherwise 


Expanding  the  spherical  harmonics  and  the  Neumann 
functions  in  terms  of  £,  n  and  <J>,  we  have: 


r=o     ^.o      $i0 


W       to 
V       V 
1        2 


m  2u   £,-m-2v.   £.-m-2v,  „,    i,-2v  -2u+v-r+2,  -2v  -s 
x  C  C   C  *      J  C  2      2  2(kR)  *    *         22 
u  v   r         s 

v+r+s     2  £  -m-2v  -r+£  -2v  -s  2  ,  , 

x  c    s  (1-5  )  *      i     2    2    {  (R/2-RA0-cZR/2) 


h*  £  (£1-2v1-2u+v-r,£2-2v2-s,C) 
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+  S/(2k)    h5        (A-1-2v1-2u+v-r-l,£2-2v2-s,C)}} 


where  the  h  -functions   are  defined   as: 
*l   £    (r.s.C)   =    HI       I   (-DU+V+kl+k2aiH,2k1+u) 
x    (&2+i, 2k2+v)Q[u^V' 5) (r+2k1+u, s+2k2+v , I ) 


and  the  Q- functions  are  defined   as 


Q£?iI'5)(P,9,C)    =    {    cosCKHS-iCAj-i^ir)    («u+v,0+6u+vf2) 


+   sin    (KRC-i(A.-&2)ir)    6    .      .(5      n-6      .,)} 
\  *\    i      z/  u+v,l     u,0      u,l 

!**()+ y*<p) 

x  cost/{tP(2i<RC-t)q}dt 

+  {sin(KR5-i(A1-A2)ir)    («u+v,0+6u+v, 2> 


+  cos(<RC-4(^-^2)tt)    5u+Vjl(<Sul-6u0) 
ricR("jH-y*cj)) 


-Ws-ty.C! 


sint/{tP(2KR?-t)q}dt 


+    {cos(kRC+4(A1+A2)tt)    <5u+v0-6u+v>2) 
+   sinURC+M^+il^TT)    6   +        } 

x 


l/{tp(2<RC-t)q}dt 
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To  conclude  this  section,  we  note  that,  as  we  did  in 
the  force  calculations,  we  have  reduced  all  the  integrals 
into  a  sum  of  one-dimensional  integrals  which  involve  the 
sine  and  cosine  integrals.  We  shall  discuss  the  evaluation  of 
these  integrals  in  Appendix  G.  Also,  we  note  that  the  actual 
calculation  is  not  as  awkward  as  the  equations  might  suggest. 
Since  the  ^-values  used  as  basis  functions  are  always  small, 
we  have  only  very  few  terms  in  all  the  summations,  as  we  noted 
in  the  force  calculations. 

5.3  Results  and  Discussion 
The  dipole  moment  calculations  have  been  carried  out 
for  several  diatomic  molecules:  lithium  hydride,  boron  hydri- 
de and  carbon  hydride.  Lithium  hydride  is  investigated  be- 
cause it  is  the  simplest  heteronuclear  diatomic  molecule, 
boron  hydride  because  it  is  more  complicated  but  still  only, 
involves  a  orbitals  in  its  ground  state  and  is  of  closed- 
shell  configuration,  and  carbon  hydride  because  it  also  in- 
volves a  n  orbital.  Also,  there  are  experimental  values  or 
some  other  theoretical  results  of  the  dipole  moments  for  these 
molecules  with  which  we  can  compare  to  our  results. 
Lithium  Hydride  LiH: 

The  ground  state  of  LiH  is  I    ,  and  the  electronic 

2   2 
configuration  is  la  2a  .  The  states  la  are  formed  from  the 

two  Is  orbitals  of  the  lithium  atom  while  the  two  outer  elec- 
trons occupying  2a  orbitals  form  the  bonding.  Thus,  the  ground 
state  is  ionic  (Li  H~). 
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(1)  In  order  to  obtain  a  dipole  moment  curve  with 
the  dipole  moment  against  the  bond  distance.  Four  non-spin- 
polarized  MSXa  calculations  have  been  carried  out  for  four 
different  separations:  R=2.75  a.u. ,  2.836  a.u.,  a. 015  a.u. 
and  3.25  a.u.   Each  of  the  calculations   is  done  with  the 
optimized  6  value   (g=RLi/R)  which  had  obtained  by  Linsen- 
meyer  (see  Linsenmeyer  (1974)),  so  that  the  muffin-tin  total 
energy  obtains  the  minimum  at  each  of  the  four  separations: 
8  =1.34/2.75  at  R=2.75  a.u.,  &=1. 39/2. 836  at  R=2.836  a.u., 
3=1.49/3.015  at  R=3.015  a.u.,  and  6=1.63/3.25  at  R=3.25  a.u. 

Also,  the  spin-polarized  a  values  are  used:  a, .=0.77159, 

Lii 

aH=0. 77627,  aii=aout=(aLi+aH)/2 ,  and  the  I   values  are  used 
up  to  3  in  the  lithium  atomic  sphere  and  the  outer  sphere 
while  £=0,1  in  the  hydrogen  atomic  sphere. The.  results  of  the 
dipole  moment  values  calculated  by  using  the  converged  MSXa 
wavefunctions  are  listed  in  table  (5.1).  The  dipole  moments 
are  plotted  against  R  in  figure  (5.5). 

The  value  of  the  dipole  moment  at  the  experimental 
equilibrium  separation  R=3.015  a.u.  is  y   (p)=  7.289  Debyes. 
Compared  with  the  experimental  value  u   =5.879  Debyes  or 
with  the  configuration  interaction  result  obtained  by  Bender 
and  Davidson  (see  Bender  and  Davidson  (1968)):  y  =6.04  Debyes 
our  result  is  about  24%  larger.  However,  if  we  look  at  the 
value  obtained  by  assuming  the  muffin-tin  charge  density  with 
equation  (3.6):  u   (p")=3.332  Debyes  which  is  about  44%  less 
than  the  experimental  value,  we  see  that  our  result  is  defi- 
nitely an  improvement. 


TABLE  5.1 


1  + 
Dipole  Moment  Components  for  LiH  (  E  ) 
at  R=2.75,  2.836,  3.015,  3.25  a.u. 
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R=2.75 


R=2.836 


K=3.015 


R=3.25 


la 

-0.000660 

-0.000547 

-0.000370 

-0.000217 

Li 

1.35972 

1.4U261 

1.49290 

1.5y927 

out 

0.U13375 

U. 011243 

0.007885 

0.004774 

la 
y 
II 

0. 

0. 

U. 

0. 

2a 

-0.540407 

-0.573004 

-0.642463 

-0.737138 

MLi 

0.020414 

0.020601 

0.21286 

0.022991 

2a 
u 
Mout 

-0.526761 

-0.5184Y9 

-0.500600 

-0.480344 

2a 
VII 

-0.22059 

-0.23700 

-0.27030 

-0.32862 

la 
^e 

1.37244 

1.41331 

1.50041 

1.60383 

2a 
^e 

-1.26734 

-1.30788 

-1.39208 

-1.52311 

total 
ye 

0.21020 

U. 21086 

0.21666 

0.16143 

yN 

2.890 

2.y48 

3.085 

3.230 

lixa(p) 

2.680 

2.737 

2.868 

3.069 

*  R  in  a.u. 
y  in  e-aQ. 
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yX0t(?O  (ea0) 


2.95 


2.85- 


2.75" 


2.65- 


2.55" 


2.8       2.9       3.0       3.1       3.2     R(a.u. ) 


Fig.  5.5  Dipole  Moment  v.s.  beperation  for  LiH  (  E  ; 
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TABLE  5.2 


Dipole  Moments  and  Their  Derivatives  for  LiH  at  R=3.015  a.u, 


p'.  ID)  Is  CD/a.u.) 


C.I.  fej.04  1.225  ,  1.31 

Exp't  5.88  1.068 

M-T  Density  3.33  2.38 

Our  Method  7  ..2a.  1.27 


a  See  Yde  and  others  (1972) 
b  See  Ebbing  (1962) 
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TABLE  5 . 3 


3  and  £-  Dependence  of  Dipole  Moment  for  LiH  at  R=3.015  a.u. 


#1 


#2 


#3 


yla 

H 

-U. 000370 

-U.0U03y9 

-U. 000370 

10 

V 
Li 

10 

y 
out 

1.49290 

1.53599 

1.4y291 

0.007885 

0.009144 

0.00/885 

ylo 
II 

0. 

0. 

U. 

y2o 
H 

-0.642463 

-0.643122 

-0.643129 

2o 
Li 

0.021286 

0.020384 

0.021401 

y2a 
out 

-U.5006U0 

-U. 486677 

-0.494427 

20 
^11 

-0.27U30 

-0.20944 

-U. 25163 

C 

1.50041 

1.54474 

1.50042 

y2° 
e 

-1.39208 

-1. 31885 

-1.36778 

total 
e 

U. 2166b 

0.45177 

0.26528 

y 

N 

3.085 
2.868 

3.285 
2.833 

3.085 
2.820 

*1:  3=1.49/3.015, 
sphere,  £=0,1 

ff2:  3=1.44/3.015, 
sphere,  £=0, 1 

#3:  3=1.49/3.015, 
sphere,  £=0,1 

y  in  unit  of  e-an 


£=0,1,2,3  in  both  the  Li  and  the  outer 

in  the  H  sphere. 

£=0,1,2,3  in  both  the  Li  and  the  outer 

in  the  H  sphere. 

£=0,1,2  in  both  the  Li  and  the  outer 

in  the  H  sphere. 
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Also,  the  value  of  the  derivative  of  the  dipole 
moment  at  the  experimental  separation  is  equal  to  1.27 
Debyes/a.u.  Compared  with  the  experimental  value   1.068 
Debyes/a.u.  and  two  other  values  by  Yde(1972)  and  Ebbing 
(1962)  (1.225  Debyes/a.u.  and  1.31  Debyes/a.u.),  our  result 
is  very  satisfactory.  Furthermore,  we  again  note  that  it  is 
an  improvement  over  the  value  (2.38  Debyes/a.u.)  which  is 
about  two  times  the  correct  value.  The  different  results  are 
listed  in  table  (5.2). 

(2)  In  order  to  check  the  sensitivity  of  the  value 
of  the  dipole  moment  to  the  way  of  partitioning  the  coor- 
dinate space,  we  repeat  the  calculation  at  R=3.015  a.u. 
but  with  a  different  value  of  B  (8=1.44/3.015=0.477).  The 
results  are  listed  in  colume  2  of  table  (5.3).  The  result  of 
the  dpole  moment  for  this  8  is  2.833  eaQ.  Comparing  it  with 
the  value  obtained  with  the  optimized  B  value  (6ODt=0.499) , 
2.868  eaQ,  we  see  the  value  of  the  dipole  moment  is  not  very 
sensitive  to  the  change  of  B,  at  least  for  this  molecule. 

(3)  In  order  to  check  the  sensitivity  of  the  value 
of  the  dipole  moment  to  the  number  of  basis  functions  used, 
we  repeat  the  calculation  at  R  =  3.015  a.u.  but  with  a 
different  set  of  bases:  I   =0,1,2  in  the  lithium  sphere 
and  the  outer  sphere;  I   -   0,1  in  the  hydrogen  sphere.  In 
both  cases  we  use  8  =  1.49/3. 015.  The  results  are  listed  in 
column  3  of  table  (5.3).  Comparing  this  result  u=2.820  eaQ 
with  the  previous  result  y=2.868  eaQ,  we  see  that  the  value 
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of  the  dipole  moment  is  quite  sensitive  to  the  basis  func- 
tions used. 
Boron  Hydride  BH: 

The  ground  state  for  BH  is   E   and  the  electronic 

2   2   2 
configuration  is  la  2a  3a  .  The  MSXa  calculation  is  carried 

out  for  this  molecule  at  the  experimental  equilibrium  sep- 
aration Re=2.336  a.u.  The  ratio  of  the  radii  of  the  spheres 
is  chosen  to  be  8=Rg/R=l .'59/2.  336  that  optimizes  the  muffin- 
tin  total  energy  at  R  .  Also,  we  use  the  spin-polarized  a 
values:  aB=0. 76206,  aH=0. 77627,  aTI=a  ut=(aB+aH)/2 •  Tne  l~ 
values  used  are  up  to  3.  The  results  for  the  dipole  moment 
calculation  are  listed  in  table  (5.4). 

The  result  obtained  is  u   (p)=-1.24  Debyes.  The  exper- 
imental value  for  the  molecule  is  not  available.  However, 
comparing  it  with  the  result  calculated  by  Bender  and  David- 
son using  natural  orbital  configuration  interaction  method: 
y=-1.47  Debyes,  we  have  quite  a  good  agreement.  Also,  the 
dipole  moment  calculated  with  the  muffin-tin  charge  density 
as  suggested  in  Chapter  III  is  uxa(p)=+1.15  Debyes  which  is 
far  away  from  both  of  the  two  results  obtained  above. 

Carbon  Hydride  CH: 

2 

The  ground  state  for  CH  is   n  at  the  equilibrium 

separation  Re=2.124  a.u.  and  the  electronic  configuration  >-■ 

2   2   2   1 

is  la  2a  3a  Itt  .  The  MSXa  result  is  obtained  with  the  follow- 
ing parameters: 

0  =  Re/R  =  1.543/2.124  which  optimizes  the  muffin- 
tin  total  energy, 
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TABLE  5.4 


Dipole  Moment  for  BH  (  Z  )  at  R=2.336  a.u. 


la 


2a 


3a 


l 

i 

^B 
i 
Pout 


'II 
i 


■0.U000015 
0.745994 
U. 0000345 
U. 
0.746023 


-0.134301 
0.084t>68 
-0.148342 
-U. 228517 
-0.42659 


-0.086499 
U. 417390 
0.6649U1 

-0.0u09y5 
U. 994996 


total   =   ^   6288:d8   e_a 


^N 


=   2.140   e-a 


0 


yxa(p)   =   -0.489  e-a0  =   -1.24  D 


Uc.iV   =   -1-47  D 


u      (p)    =   +1.15   D 


See  Bender   ana  Davidson    (1969) 


TABLE  5.5 


Dipole  Moment  for  CH  (  II)  at  R=2.124  a.u, 
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la 


2a 


3a 


■ITT" 


-0.0o0000l4  -0.042175  -U.0b24l9  -0.000543 

0.581002  0.174046  0.3k»7942  0. 259628 

0.00000380  0.025876  0.288481  0.171516 

0.  -0.115285  -0.089127  -0.033806 

0.581005  O.0424S3  0.534878  0.396796 


l 

W 

y1 
ii 
i 


V     =  2.713488  e-aQ 

yN  =  1.943  e-aQ 

yx  (p)  =  -0.7705  e-aQ  =  -l.y58  D 


^C.I 


=  -1.472  D 


y   (p)  =  -0.04  D 
^exp't  --l-«6  +  0.06  D 


See  Bender  and  Davidson(ly69) . 
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a~  =  0.75331;     a  =0.77627;     aTT  =  a   ,.=0. 76479  ; 
C  h  II    out 

£  =  0,1,2,3;     £=  0,1    ;     I  =0,1,2,3     ; 

C  H  out 

for  the  la,  2a,  3a.  states 
£c  =  1,2,3   ;    %=  1,2    ;    £out  =  1,2,3       ; 
for  the  Itt  state. 
The  results  of  the  dipole  moment  calculation  are  listed  in 

table  (5.5) 

a, 
The  result  obtained  is  uxa(p)=-1.95«  Debyes.  Compar- 
ing it  with  the  experimental  value :u   =-1.46+0.06  Debyes, 
and  also  with  the  theoretical  value  obtained  by  Benaer  and 
Davidson  using  natural  orbital  configuration  interaction, 
urT=-1.472  Debyes,  we  note  the  value  that  we  got  is  not  too 
different  from  these  values.  Indeed,  it  is  a  big  improvement 
over  the  value  obtained  by  using  the  muffin-tin  charge  den- 
sity: yXa(p)=-0.04  Debye. 


CHAPTER  VI 
CONCLUSIONS 

The  main  purpose  of  this  dissertation  was  to  investi- 
gate the  quality  of  the  MSXa  wavef unctions.  To  achieve  this, 
we  performed  the  force  and  the  dipole  moment  calculations  by 
use  of  the  charge  density  generated  from  the  MSXa  wavef unctions ; 
instead  of  using  the  muffin-tin  part  of  the  density  function. 

With  the  results  obtained  for  the  dipole  moment  of 
LiH,  BH  and  CH,  we  see  that  the  use  of  the  MSXa  wavef unctions 
in  such  a  dipole  moment  calculation  is  quite  adequate  even 
though  it  is  well  known  that  the  dipole  moment  is  very  sensi- 
tive to  the  approximate  wavefunctions  used.  The  reason  that 
we  may  expect,  and  we  do  have,  adequate  results  for  the  calcu- 
lation of  the  dipole  moment  with  the  MSXa  wavefunctions  is 
as  follows.  .  In  the  MSXa  method,  the  constant  potential  and 
charge  distribution  in  the  intersphere  region  is  the  most 
severe  approximation,  and  thus  the  MSXa  wavefunction  differs 
most  from  the  exact  Xa  wavefunction  in  this  region.  However, 
the  amount  of  the  total  charge  in  this  region  is  much  smaller 
than  that  in  the  other  regions, and  furthermore,  if  we  take 
the  center  of  the  outer  sphere  as  the  origin  of  the  whole 
system,  we  have  quite  a  large  cancellation  of  the  errors  in- 
evolved  in  the  contribution  by  the  electronic  charge  cloud 
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Ill 

in  the  intersphere  region.  Therefore,  the  net  result  of  the 
dipole  moment  calculation  is  quite  goad  with  such  MSXa 
wavef unctions. 

In  the  treatment  of  the  heteronuclear  diatomic 
molecules,  the  calculation  is  more  complicated  because  here 
we  have  an  arbitrary  choice  of  the  parameters:  the  ratio  of 
the  radii  of  the  spheres,  and  the  a  values  in  the  intersphere 
region  and  in  the  outer  sphere.  Though  we  have  not  broadly- 
investigated  their  effects  on  the  results,  we  feel  that  they 
do  effect  the  result,  as  we  have  noticed  in  the  treatment  of 
LiH. 

However,  in  the  treatment  of  homonuclear  diatomic 

molecules,  we  do  not  have  such  a  problem  of  the  choice  of  the 

a  values  and  the  3  values.  The  two  atomic  spheres  have  to  be 

the  same  size  and  the  a  value  has  to  be  a  constant  over  the 

whole  space.  Thus,  in  the  Hellmann-Feynman  force  calculations 

for  H2  and  N2,  there  are  not  any  adjustable  parameters  except 

the  a  value.  Since  the  Hellmann-Feynman  Theorem  is  satisfied 

by  the  exact  Xa  wavef unctions  no  matter  what  a  value  is  used, 

the  only  approximation  that  causes  the  electrostatic  force 

to  differ  from  the  quantity  — jLXct >   is  the  muffin-tin  approx- 

imation,  through  which  the  MSXct  wavefunctions  u.^  are  obtained. 

From  the  results  for  H2  (f igure(4.4) ) ,  we  notice  that  there 

is  such  a  difference  between  the  two  quantities  Fxa(p)  and 

ii      . 
_8<EXa(p)>   Thus,  at  least  for  the  case  of  H9,there  is  a, 

3Xp 
difference  between  the  MSXa  wavefunctions  and  the  exact  Xa 

wavefunctions. 
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As  stated,  another  purpose  of  this  dissertation  is 
to  suggest  a  method  for  the  evaluation  of  the  complicated 
three-dimensional  integrals  over  the  intersphere 
region.  Actually,  this  method  of  using  the  Prolate  Spheroidal 
Coordinates  should  not  be  limited  to  the  force  and  dipole 
moment  calculations.  It  can  be  employed  in  the  evaluation 
of  molecular  properties    that  can  be  represented  by  some 
one-electron  operators.  Thus,  for  example,  the  transition 
amplitudes  <u-|r|u>  between  any  two  molecular  states  should 
be  able  to  be  evaluated  by  this  method.  Of  course,  the  result 
would  depend  on  whether  the  MSXa  wavef unctions  are  accurate 
enough  for  the  evaluation  of  such  a  sensitive  quantity.  The 
merit  of  this  method  of  calculation  is  that  all  the  integrals 
are  reduced  into  some  one  dimensional  integrals  that  can  be 
handled  in  an  efficient  way. 

Although  in  this  dissertation  we  have  only  treated 
diatomic  molecules,  we  can  also  apply  our  method  to  more 
complicated  systems.  To  illustrate  this,  we  give  a  brief 
discussion  on  the  calculation  of  a  quantity  F  for  a  triatomic 
molecule  ABC. 

As  before,  we  partition  the  space  into  the  outer 
region,  the  intersphere  region  and  the  atomic  region( in  this 
case,  it  includes  the  three  spheres  A,  B,  C)(f igure(6. 1) ) . 
The  difficult  part  of  the  calculation  is  in  the  evaluation 
of  the  following  integrals  over  the  intersphere  region  II: 

F1  =  f  u2(r)F(r)d3r 


113 


out 


Fig.  6.1  Schematic  Representation  of  a  Triatomic 
Molecule  in  tne  MSXa  Method 
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*\»  -*  '  th 

where  u.(r)  is  the  converged  MSXa  wavefunction  for  the  i 

molecular  state. 

As  before;  the  evaluation  of  F1  can  be  reduced  to 

the  evaluations  of  integrals  like: 


F 
AB 


1   =   f   n   (<r  )Ymi(r"  )n   (icr  )Ym2(?  )F(?)d3r        (6.1) 
AB    J    £i    A  l1     -AJ    9.2         B   £2   B  ■  K         J 

We  can  evaluate  F^g  by  changing  the  region  of  inte- 
gration to  the  intersphere  region  plus  the  sphere  C,  and 
then  substract  from  it  the  contribution  from  the  atomic 
sphere  C: 


'In-     I 

'JI  +  C 

-  I 


^/KrA>V.rA>"*2^B>Yi;<rB)F<r>d  r 


^,(>;rA)<I(?A)\2(KrB)Y":(fB,F<f)d3r    <6-2) 


To  evaluate  the  first  term,  we  note  that  since  the 
integrand  has  poles  only  at  centers  A  and  B,  and  both  of  them 
are  excluded  from  the  region  of  integration,  there  is  no 
problem  of  divergence.  So  this  term  looks  like  the  previous 
integrals  that  we  encountered  in  the  diatomic  case,  and  we 

want  to  apply  our  method  to  the  evaluation  of  this  term. 

m  -»■ 
However  ,    generally  the  spherical  harmonics  Y^(r.)  and 

Y^(r  )  are  not  defined  relative  to  the  internuclear  axis  R^B 

as  they  are  in  the  case  of  a  diatomic  molecule.  Thus  we  first 

rotate  the  coordinates  of  both  A  and  B  in  such  a  way  that  the 

new  z-axes  of  both  atoms  are  oriented  along  the  internuclear 
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ZA 
A- 


XA 
? 


Fig.  b.2  Rotation  of  Atomic  Coordinates  to  Align 
with  the  Internuclear  Axis. 


lib 

axis  R«B  and  both  the  x-axes  and  y-axes  be  coplanar  with  RAj 
(f igure(6.2) ) .With  this  rotation,  the  complex  spherical  har- 
monics are  transformed  as: 


i 

I 

,'=-1 


v"(vv  -i*?'<°k-*k>  dL'(v8a-v 


where  D   ,  are  the  rotational  matrix  elements  and  a»,  &a,    Ya 
mm'  aaa 

are  the  Euler  angles  of  the  new  coordinates  with  respect  to 
the  old  coordinates  centered  at  A.  Since  the  real  spherical 


harmonics  Y„  can  be  written  as 


\    =  ^l    +    {~1}    If  I    )/2  f°r  m>° 

the  transformed  real  spherical  harmonics  are: 

Ym(e  ,d>  )  =  T    Hm'(6'  ,<J>'  )  JQ.  l      (u>     )         (6.3) 
V  A'%     t-.  ,Q  I    V  A,HA  M   mm'   A.AB' 

o  j,  m  % 

where  Z^M^   =  Dnmi'<aA'BA^A>- +  ("^  DEm'(aA'BA^A) 

and  similarly  for  Y^CBg,  <J>B) .  Quite  similarly  is  for  m<0. 
Thus  the  first  term  in  equation  (6.2)  is  equal  to: 

2_,   ^/m,m'   A,AB   *mm  '   B,AB 

and  now  we  can  apply  our  method  to  the  evaluation  of  these 
two-center  integrals.  Before  that,  we  have  to  expand  F(r) 
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about  A : 

CO        <1 

F(r)  -LI,      F  3(r'^  3(r') 

4^=6   vn^-Jtj    3        3 

Then,  we  wind  up  with  the  evaluation  of  the  integrals: 

JIJ  +  C     1  1  2  2  3  3 

By  use  of  the  Prolate  Spheroidal  Coordinates  with  the 
foci  at  A  and  B,  we  can  reduce  the  integral  into  a  two-dimen- 
sional integral  through  the  evaluation  of  the  <j>  coordinate: 

>=  2*  Vm^.O  j[i<rt  n£l(<rA)P^(coseA)n,2C<rB)P^(cosy-) 

-mJ-ml  -m'-m'  3 

x   Fu    *      2(r    )P„    *      2(cos9*)    d  r/d<j> 
J63  A      J63  a 

Since  PA1(cos9')P    2(cos9')P~    r~    2(cosU') 

£j  A      £2  B      £3  A 

m.'    m!    -m;-ml  m'  m'  -m'-m' 

=  N,1^    2N„    x      2(.sin9')    l(sin8J.)    2(sin9')      a      2 


r  x:     £  s: 


01  0) 


X     (C0S9')      1113  12  3      CcosQ-)      2         2  2 


and  (sine*)    1(smQ')    2(sin9')      -1      2 


. , mi ,  .    -ml 

=    (sin9')    2(sm9')      2 


118 


{(S2-i)(i-n2)>imV(5-n)m;   *   wf-lHl-S)}-*"''/^)** 


m' 
m+n)/U+n)}    2 


ana  similarly  for  the  cosine  factors.  Thus,  we  can  apply  our 

method  of  integration  as  before,  except  that  the  limits  of 

integration  are  more  complicated. 

For  the  second  term  in  equation  (6.2),  we  can  expand 

n   (<r  )Y°1i(r  )  and  n   Or  )Ym2(r  )  about  center  C.  Since 
A/   A   £x   A       *2    B   £2   B 

inside  sphere  G,  we  have  rp<R.p  and  r-,<  RRr,   so  the  expan- 
sion is : 


n(i<r  JY°(r  )  =  4tt  Y>     X   i^'^'Krm'ilmjrni") 

£  a  i    a  itX'   r,w 

X'V(KRAC)Yr»(5AC)J^lKrC)Yr(?C) 


m  ->-  -*■  ■ 

and  similarly  for  n„  (KrD) Y.cr-) .  Also,  we  expand  F(r)  about 

center  C: 


•*■  r^    m      m  -*■ 

F(r)  =   £  FJl(ru)YJl(rc) 


Then  we   have   the   integrals: 


Jv  „ni '    -»■  .         .  „m  '->•«_■     ,    -     „m    ,  -»■    .  ,  3 

c    J£^rc)Yr(rcJJr(Krc)Yi'trc)F£3(rc,Yi3(rc)d  r 

rRt  2 

=    ia'm';£'m';£3m3;  j£  -  (<rc}  J  £  -  UrC)F£  3  (rc}I"c   dr< 
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which  involves  only  one-dimensional  integtals. 

The  formalism  given  above  is  quite  general  for  any 
function  F(r).  In  particular,  for  the  dipole  moment  F(r)=z 
(or  x,  y),  and  the  expansion  in  spherical  harmonics  then  has 

terms  only  up  to  1=1.    Similarly  for  the  Hellmann-Feynman 

'  -*■        2 
force  function  F(r)=cos9/r  which  corresponds  to  the  spherical 

harmonic  term  1=1. 

In  conclusion,  the  use  of  the  charge  density  generated 
from  the  MSXa  wavefunctions  does  improve  the  results  of  dipole 
moment  and  force  calculations  over  the  use  of  the  muffin-tin 
charge  density.  From  the  results  of  these  two  kinds  of  cal- 
culations, we  see  that  the  MSXa  wavefunctions  are  quite 
adequate  in  such  molecular  property  calculations.  Any  further 
improvements  should  come  from  the  correction  to  wavefunctions 
that  is  a  solution  to  the  muffin-tin  approximate  Xa  one-elec- 
tron Schrodinger  equation  (2. 6) .  In  fact,  such  an  approach  is 
being  carried  out  by  Connolly,  Williams  and  others. 


APPENDIX  A 


DERIVATION  OF  THE  MUFFIN-TIN  Xa  ONE-ELECTRON  EQUATIONS 


In  Section  2.3,  we  mentioned  that  the  necessary  and 


sufficient  condition  for  the  muffin-tin  total  energy  <E   ("p)> 
be  a  ^minimum  is  that  the  spin-orbitals  u.  satisfy  the  muffin- 
tin  one-electron  equations  (2.6).  We  are  going  to  show  it  in 
this  appendix. 

Starting  from  the  non-muffin  total  energy  expression: 


where  u-  are  the  spin-orbitals, 

U„„  is  the  nuclear-nuclear  interaction  energy, 

and    Uj  is  the  sum  of  the  nuclear-electron  coulomb  energy 

U   ,  the  electron-electron  coulomb  energy  U„„  and  the  ex- 
ne  ee 

change-correlation  energy  U   : 


IL,  =  U   +  U   +  U 
T    ne    ee    ex 

=  J  P(?l)Vne(?l)d3rl  +  I  I  P(?1)Vc(P.?1)d3r1 
+  Ca  J  {^(r^)4/3  +  p^(?1)4/3}d3r1  (A.l) 

120 


121 

(-21    ,,  ,..  ..  , 
p"  '  1   p1 


where   -V   (?x)  =  £  (-2Z  )/  l^-R  | 

pal  "  ^ 


9  1/3 


It  can  be  shown  that  for  any  two  functions  f  and  g,  with  the 
corresponding  f  and  g,  we  have  the  following  identity: 


J  f(r")g(r)  d3r ■=  J"  f(r)g(r)  d  r 


Hence,  we  have 


Une=  1  p(?l)Vne(?l)d3rl=   J  p<?l)Vne(?l)dV       (A'3) 
Uee  =  2  ]VP(r1)Vc(p,r1)d3r  =|  J  p(r1)Vc(p ,r1)d3r1   (A.4) 


Now,  let  us  look  at  the  variation  of  <E  ("£)>  with  respect 
to  u.  which  is  of  spin  s  (either  up  or  down).  By  using  the 
functional  derivative  notation  (see  Volterra  (1930)),  we  have 


6  <E   (-)>  =  f  ^r-^-    <ExJp)>  6^i(r)  d  r  +  (complex  con  j  .  ) 
xa  p     J  6u  (r)    Aa 


Here  u^(r)  is  an  arbitrary  variation  of  the  function  u.(r) 
over  its  whole  range  of  definition.  There  is  no  contribution 
from  the  nuclear-nuclear  interaction  term  since  it  is  not  a 
functional  of  the  spin-orbitals .  As  to  the  kinetic  energy 
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part,  we  have; 


6u 


K.E 


^*(?y 


n.    J  uj(r1)(-V1)uj(r1)d,:>r16u.(r)d' 


+  (complex  conjugate) 


Ji 


ni6(r-r1)(-V1)ui(r1)6ui(r)  d  r^d  r2+  (c.c) 


Integrating  over  r1 ,  we  have 


5U 


K.E 


t  O/*  ->•     o<\,      ->     ? 

=  J  n.5u.(r){-?u.(r)}  d  r  +  (c.c.) 


(A. 5) 


The  variation  of  the  term  U„  is  more  complicated 
because  it  depends  on  the  spin-orbitals  only  through  the 
muffin-tin  charge  densities  p.  and  p.  and  we  have  to  use 

T  + 

the  chain  rule  of  differention  generalized  to  functional 
derivatives  to  find  U  : 


5Um  = 

T 


6UT     ^*  _>.    3 
-- — —  6u  (r)  d  r  +  (c.c.) 

<$u.(r)     x 

l 


=  Y   f  d3r  sT.ih    f  dVJgr__  *Ps'(?'> 


6ps,(r')  5ui(r) 


where   s  '  =  +  ,  4- 


+  (c.c.) 

(A. 6) 


Since 


Ps.<r')  = 


1_ 

4tt 


n .u  .(r ' )u . (r ' )dft'      if  r'  in  sphere 
J  J     J 


J.Sj**' 


L  xj —  I   >   n  .u  .  (r '  )u  .(r  '  )d  r'     if  r'  in  region  II 
Vl1  J«j^'J  J     J 


1^3 


Hence  rj7  n.{   U(r-r')u.(r')d'n'}    <5  ,       if   r'    in    3 

6Ps'(?')  Si'S 

V^Y  ni(j6(r-r")u    (r")d   r"}5c 


<Sui(r) 


(A. 7) 


,       if   r'    in    II 


Return  to  (A. 6),  let  us  consider  the  variation  of  the  term 
Une'  which  is  included  in  UT,  with  respect  to  P"  ,  .  From  equa- 
tion (A. 3),  we  have: 


<SU 


ne 


6ps,(?') 


5ps.(r«)  J 


P(rl)Vne(rl)d  rl 


«P„.<* 


^  I   J   Ps(?l)Vne(?l)d  rl 


=  J   f  gPs'(ri)  y   (?  )d3r 

r  J  5PS.(?')  ne  x     x 

-  iK's  6(?i-?,>  w  d3ri 


=  V  (r») 
ne 


(A. 8) 


Secondly,  let  us  consider  the  variation  of  the  elec- 
tron-electron coulomb  energy  term  U   with  respect  to  the 
variation  of  p~gf.  From  equations(A.  1)  and  (A. 2),  we  have: 

SUee 


6p_,(r«)     5p  ,(?•) 


_  1  [  p(r  )V  (p,r  )  d3r 
<  )    2     J  1   c    1     1 


^7?~  2    ?  I  W    2PS<?1>  V(?2>/rl2   «  V 


"      ^    £  j  ^^    's(V     j    <"' 


"2r12      ps.(f') 
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d  r2  Ps„(r2) 


3    1   6Ps(ri) 


J  A, 


ri2  «Ps'(r') 


s  ^° 


12 


Xl[d~r2ps„(r2)  jd^^-  iUj-f')  «ss, 


S   S" 


•  3   JL 

'12 


=  I   d  r   p  (r  )/  r  -r» |  +    L        &   r  p  „(r  )/|r'-r 
j  J     1   s   1     1        s»  J  •   2   s"   2 

,  Jd  ri-2ps(ri)/ivr,i 

J  d3r]L  2  .p(?1)/|?1-?'  |  =  Vc(p,r') 


(A. 9) 


Finally,  let  us  consider  the  exchange-correlation 
energy  term  U   .  From  equation  (A.l),  we  have: 

6X 


6% 


.fips.(r')    5ps,(r') 


f   <--/-►  n4/3    3 


's   1 


t-  f       —   -*■   1/3  „  -»■   -»■ 
-  Ca  I  "I  (4/3)  Ps(ri)     6(rrr' 


-   -»■   1/3 
=  (4/3)  Ca  Pgf(r')  ' 


)  6    dr. 
ss'     1 


(A. 10) 


Combining  equations  (A. 8),  (A. 9),  (A. 10)  and  substituting  into 
equation  (A. 6) : 


6U„ 


=  X    j    d3r      u.(r)     jd3r'    (Vne(r"'  )+Vc(  p',?'  )+|cciPs  -  (r1  )1/3> 
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6ps,(r')' 

x   -^rr:   +  <c-c-)  (A.ii) 

<5u. (r) 
l 

_ .  —   -*■  o<*  ->- 

Since  we  have  two  different  forms  for  6p  ^t'  )/Su..(r)   depend- 
ing on  the  location  of  r' ,  as  expressed  in  in  equation  (A. 7), 
we  separate  the  integral  over  r'  into  two  parts,  one  over 
the  atomic  regions  ( 6=0, 1 ,...  N)  and  the  other  over  the  in- 
tersphere  region  II: 

6U  =  Z   d  r  6UjL(r){  L        dV  — V4"  {Vne(r«  >+Vc<P  • rR> 
s'  J  (J"-°  J?        <5uj_(r)    ne   6    c     6 

3  6P_,(r')   ._ 

i   d  r>  — -§- {  V   __ 

in  ^*/-*-\             ne,II 

11  6u.(r) 


4  „     -  ,+    ,1/3,     f    3    <5P„,(r'  ) 
+  fca  Ps.(rs)    >■+   ]   d  r-  — S_ {v 


+  Vc(p,?')  +  |  Ca  PII>S,}}  +   (c.c.) 


By  equation  (A. 7),  we  have 


OTT 


=  X  j  d3r  u*(r){£  [d3r'  {  l/(4ir)  ni« 

x    SCr'-r^u^r''^}  ^ne(rg)+Vc(P,?3)  +  (4/3)CotPsl(?')1/3} 

+  Jn  d3p'{V^  ni  6Sis>  J/t^'V^^"1^,!! 

+  V  (p,?')+(4/3)C  pTT   ,1/3}}   +   (c.c.) 
c  a  II , s ' 


Integrating  over  r,  we  have 
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6Un 


N  f 

X    J    d3r^    (1/4tt)    n.    {   j    6u*(r  •  )u\(r  •  )d$2£   } 


pro    Vp 


*    {Vne(r6)+Vc(^?e)  +  (4/3)CaPSi(r3)1/3} 


+     [     dV    r^-  n,        j     6ut(r")ui(r,)d3r" 

x    ^ne>II+Vc(p,?')+(4/3)CaFII    s1/3}+      (c,c.) 

'    i 

=    I   [r82drBdfl6  "i/<**><   Vne(r')+Vc(p,;')+  |caPs/rB)1/3} 

(3  =  o    ^  i 

f  ^*    -»■      ^      -»•  I         3  %*    -»■      a,      ->■         1 

x{  dfi'    fiu-d-JUifr')}   +  d  r"    6u.(r")u.(r")-^-  n. 

J  P  lfc!-LP  ■'II  1  1  VH1 

x     j„   dV{  Vne,II+VP'?'>+  'KpiifSi1/3>      +      (c.c.)-- 
=     X    V(4.)  j   r^    {   J 

+    (4/3)C   p      (rl)1/3    }x{  f    dfi>6u.(r')u.(r')} 

53  j_        p  *■  p        1        p        1        p 


dB8     Vn,(ri>+ve<P.'B> 


1/  -►  1<  -♦■ 


{3  -v*     ->        ^        -►  1  f 


<Vne,II+Vc(P.'r,) 


+    (4/3)C    p 


1/3 


aII,si 


}  +      (c.c.) 


(A. 12) 


Since  the  muffin-tin  average  of  a  function  f  is  defined  as: 


f(r)  = 


1/(4tt)J  dft  f(r)  if  r  in  the  atomic  regions 

f   3    ■*■  -*■ 

1/Vjj   d  r  f(r)  if  r  in  the  intersphere  region 

Jjr 
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So,  equation  (A. 12)  is  reduced  to 


<SUm  =  y        n./(4ir)  [  r '  dr '  {4ttV   (r')+4TTV  (p,r') 


T    _. 


+  47r(4/3)CaPs    (r^)}{     J     d^'    u.(r')    u*(?^)} 

+  VVII    h*3*"    &1&'A&'>    {Vll\e,II   +VIIVc,II(^ 
+  4tt(4/3)C    pTT         }  +  (c.c.) 

=       X      n  d3r^    «3I*(?.){   Vne(r')+Vc(p,pj> 

+    (4/3)CaPSi(r^)1/3}   ^(^)    +   n±  j^dVfi^'H    VneII 

+  Vc,II   +   (4/3)CaPlI,si1/^}   ^(r")   +   (c.c.) 

Therefore,    <5UT  =  \  d  r  n.  6u. (r){  V  (r)  +V  (p",r) 
i    )  i   i      ne      c 


+  (4/3)CaPg_(r)1/3}  u.(r")  +  (c.c.)    (A. 13) 


Thus,  combining  with  the  kinetic  energy  part  in  equation 
(A. 5),  we  have: 

6<E   (p)>  =  J  d3r   n.  u,(r)  {-V2  +V   (r)+V  (p,?) 
xa  y  J  x      1  ne     c 


4  —   -»•  1/3  ^  -»■ 
+  |Caps.(r)    }  ui(r)   +   (c.c.)  (A. 14) 
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Since  we  have  the  restraints  on  u.  imposed  by  the 
ortho normality  conditions,  we  have  to  introduce  the  Lagrange 
multipliers   A.  .   and  demand  that 

6  {<Ev    (p)>  +   XX  n*'n*  Ai  ,   [  u?(r)u.(r)  d3r}  =   0 
xct  M       •,  j   1   j   iJ  J   1     j 

XX 

where  we  introduce  a  factor  n?n^  with  each  A.  .  for  the  con- 
venience of  the  derivation. 

Since   <5  {  £  X  n  "   *     u.(r)u.(r)  d  r} 
i   j   i  J   ij  J   1    J 

=   4.  ni  nj  Aij  J  <Sui(r)uj(r)ddr   .+    (c.C.) 
j 

therefor,  together  with  equation  (A. 14),  we  have : 

j  d3r  6u*(?)  {n.{  -V2  +  Vne(?)  +  V^,?) .+  \   Caps.(?)1/3} 

x  u  (r)  +  nf  "T  n*  A.  .  u.(r)}  ■  +   (c.c.)   =  0     (A. 15) 

1     ■   i  y   J   iJ   J 

If  the  first  part  of  equation  (A. 15)  is  equal  to  zero, 

so  will  be  its  complex  conjugate.  Thus  we  can  drop  the  com- 

o,*  ->■  ■ " 
plex  conjugate  term.  Also,  since  5u.(r)  is  arbitrary,  there- 


fore, the  necessary  and  sufficient  condition  for  <E   (n")>  to  be 

xa  M 

stable  is  that  u.  satisfy: 

{-V2+Vne(;)+Vc(p,r)+  gC^C?)1/3^^?)  =  -  2  Aij4SJ(?) 

(A. 16) 
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Since  the  \^i   matrix  is  hermitian,  it  may  be  diago- 
nalized  by  some  unitary  transformation  C: 

»!  »i  -  X  cij  n3  uj 

J 

or,  in  matrix  form:   ti  u'  =  C  n  u 

It  can  be  easily  shown  that  the  density  formed  from 

u'  is  invariant  for  such  unitary  transformation.  Hence,  the 
1 

operator  H(p)  on  the  LHS  of  equation  (A. 16)  remains  as 

j. 
the  same  form  after  the  transformation,  i.e.,  CHt(p)C  =H(p'). 

Writing  equation  (A. 16)  in  matrix  form: 


H(p)  ffi  v.  =   -   )\   n~  <u 


Multiplying  both  sides  by  C  and  insert ing  CtC(whichis  equal 

i 
to  identity  because  C  is  unitary)  in  between  H  and  m2  on  the 

LHS  and  in  between  \\  and  n2  on  the  RHS ,  we  have: 


t    A  %  t    i  ^ 

C  B(p)  C  C  <n2  u  =  -  C  X  C  Cm2*! 


which  implies    H(p')  ti  u'  =  -  V  u  u' 


a. 
Now,  since  V  is  diagonal,  therefore  A!  .=  5.  .  e!.  Thus 

O  _     -»■   _   Li.   ->     4    —     -*■  1  /  3 ,   ^   ,  -*■ „     ^  .  ^  .  ,  "*"  >. 

{"V  +V   (r)+V  (p',r)+  ^C  p'  (r)    }  u'(r)  =  £;u  (r) 

lie         C  O  Ut  O  A  J-  j_   X 

which  is  the  muffin-tin  Xa  one-electron  equations. 


APPENDIX  B 

THE  MULTIPLE  SCATTERING  METHOD 

To  solve  the  muffin-tin  approximate  Xa  one-electron 
equations,  we  can  employ  the  multiple  scattering  method  which 
is  commonly  used  in  scattering  problems.  Starting  with  the 
equation: 


{-V^  +  V   (p,r)}  u.(r)  =  e.u.(r)  (B.l) 

xa        1       ii 


we  have,  in  the  atomic  regions: 


(-'2  +  vx0(p'V}  w =  wv 


The  solutions  inside  the  sphere  $   can  be  expanded  in  a  par- 
tial wave  expansion: 

VV  "  X  cj£  VVV  Y-(?6)  (B.2, 

where  we  have  used  the  real  spherical  harmonics  and  the  ra- 
dial functions  R„  have  to  satisfy  the  following  differential 
equations: 
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{"  -\  T~   (r!   rtT")    +   MA+D/rj!   +Vv    (P.r.)}    Rff(rfl,e    ) 
rR2   drft        B   drR  B        xa  B  I      3 


=   e      R   (r    ,e    ) 
i      £      8      i 


In  the  intersphere  region  II,  the  potential  is  a 
constant  -V-,.  Thus  equation  (B.l)  is  reduced  to: 


I  -V  +  (V   -  e  )}  u  (r)  =  0 
II     i    i 


which  is  just  the  equation  for  a  free  electron.  Hence  we  may 
express  u-  in  this  region  as  a  multicenter  partial  wave  ex- 
pansion: 


r  T   T  Ai6n  (<r  )Y?C?')  +  I  A 


a.  -»■ 
u.(r)  = 


10  •  ,  \vm,-+      X 


(B.3) 


>   V  iB        m  -*     V  iO 


JmV^VVV 


if .  e  <  V 


II 


where    <  =  ( |  V  -e . | ) 

n«,  jg    are  respectively  the  ordinary  spherical  Neu- 
mann and  Bessel  functions:  k^  and  i^  are  respectively  the 
modified  spherical  Hankel  functions  and  Bessel  functions. 

With  the  requirements  that  the  spin-orbitals  u.  and 
their  first  derivatives  have  to  be  continuous  on  the  surfaces 
of  the  different  spheres,  the  C-coef f icients  and  A-coef f icents 
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can  be  determined.  The  method  is  to  express  each  of  the  par- 
tial waves,  n^Cicr  )Y^(r  ),  Jjl(Kr0)Yii(?o)  etc'  of  different 
centers,  in  terms  of  partial  waves  all  centered  at  one  sphere 

'V. 

a,  and  then  join  smoothly  with  the  u.  as  expressed  in  equa- 
tion (B.2)  inside  the  sphere  a,  and  similarly,  for  the  deriv- 
atives. The  result  is  that  the  C  and  A-coef f icients  are  re- 
lated by  the  following  equations: 


ia 
L£m 


L10 


Kb  [j0(Kb  ),R„(b  ,e.)]  c] 

a  L  £   a    &   a   i  J   * 


ia 
m 


ei'>  VII 


L  (-1)£+1  Kb*  [i,«ba),R£(ba,^.)]c^,  \    <  vn 


^olWV'^^Vl  c 


iO 

im 


i+1     ,  2 


iO 


(B.4) 


%  — 


L  (-i)"   <b0  [a^bo^p.kj^b,,)]  cim    ei   <  Vu 


where  b  are  the  radius  of  sphere  a, 
a  f> 

[u,vj  are  the  Wronskian  defined  as  u  -r—  -v  -r—  . 
Also,  we  obtain  a  set  of  equations  for  the  A-coef f icients: 


%   —1  i 


=   Y  ^    G        (R  „,e  )An  ,„, 


(B.5) 


where  R  0   are  the  distance  between  the  centers  of  sphere  a 
and  sphere  8;  and  the  t  and  G-coef f iceints  are  defined  as 
below: 
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t£(E)    = 


- 1 VKVA<Ve)J/K<KV«R*(Ve)1      e  >  vn 


Ve)  = 


•-  [^(■cb0).R£(b0,OJ/[iJl(xb0),RJl(b()fe)J  e    <   Vn 

(B.6) 
r(l-6    0)4TTii_J''  y,    i_Ji"    I(£"nr;£m;£'m') 


£m,£'m'      aP 


x    nt..(«Bo8)Y']l.,(BoB) 


e   >  V 


II 


(1-6    0)47r(-l)£+£'   V     Kl'mMmiiV) 
a  p  t-r** 

i' » 


m" 


4TTi)l~S''   T*      i~5,"r(£"m";£m;£,m' ) 
1  *» 


e    <   V 


II 


ao  , 

Giim,£'m'(RaO'£)    1 


m"   -»■ 
x  Jr.(KRa0)Yr,(Ra0) 


e   >  V 


II 


4tt(-1)£+£? '  "2     I(£"m";£m;£'nT  ) 


l    x    ifl„(icRrtn)Y'"    (R   n) 


m" 


e    <   V 


II 
(B.7) 


where    I(£"m" ; £m; £'m' )    are   the  Gaunt   coefficients   defined   for 
the   real    spherical   harmonics   as: 
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j    tf. 


I(£"m'  ;£m;£'m')  =   |   Y^  „(r)Y™(r  )Y* ,  (r )  dQ 


From  equation  (B.5),  we  see  that  we  have  nontrivial 
solutions  for  the  A-coef f icients  and  hence  for  the  C-coeffi- 
cients  only  when 


det  It^Ce  )  X6   6    5    -  G?6  „   , (R  ■  ,e,)|  =  0       (B.8) 
1  V    i    ct8  HV    mm'     lm,l   m'v  aB'  x  ' 


Thus  we  can  find  the  eigenvalues  e.  by  finding  the  zeroes  of 
the  determinant  (B.8),  and  from  this  we  can  find  the  C  and 

A-coef f icients ,  and  hence  the  wavefunction  u. ,  by  equations 

a 
(B.5)  and  (B.4).  Here  we  notice  that  the  functions  t«(e) 

,%  -      .* 
depend  on  the  value  of  K=|e-Vjj|"  and  the  value  of  the  radial 

function  R»  and  its  derivative  on  the  boundary  of  the  sphere 

and  therefore  implicitly  on  the  potential. The  G-coeff icients 

also  depend  on  e  and  the  potential  through  the  quantity  k. 

Hence  the  determinant  in  equation  (B.8)  is  nonlinear  in  e. 

The  way  we  determine  the  eigenvalues  e.  is  to  compute  the 

determinant  over  a  range  of  e  that  bracketed  the  guessed 

'V, 

eigenvalues  e^.    The  value  of  the  determinant  changes  sign 
on  the  two  sides  of  the  guessed  e-.  Then  we  inversely  inter- 
polate the  values  to  find  the  zero,  and  hence  the  eigenvalue 
e^,  of  the  determinant. 


APPENDIX  C 


HELLMANN-FEYNMAN  THEOREM  AND  THE  Xa  WAVEFUNCTIONS* 


In  Chapter  II  and  Chapter  III,  we  mentioned  that  the 
exact  Xa  wavef unctions,  for  the  case  where  a  is  constant  ev- 
erywhere, rigorously  obey   the  Hellmann-Feynman  theorem.  We 
present  a  proof  in  this  appendix. 

Starting  with  the  expression  for  the  exact  Xa  total 
energy  in  equation  (2.1): 

<Exa>  =  T   ni  \^1n-V1)ni^1yd3v1 

J     fr      I'l'RpI  1         X     2   ■>>  r12  1 

+  I    [  "^\a^l)d\   +  \    j    <V?l>Uxa+(?l)d3rl 


+  5     ^    2ZpZq/Hpq  (CD 

where  the  energy  is  in  Eydbergs  and  the  distances  are  in 
atomic  units.  Also, 


U™  «<r)  =  -9a  {  3/(4tt)  p  (?)  1/3  }       s  =  t  and  4- 


*  The  main  reference  for  this  appendix  is  J.C.  Slater  (1972b) 
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'(r)  =  P  (?)  +  P   (?)  =  y  n.u.(r)u.(?)  (C.2) 

t  +         *r"   i  J.      J- 


Now  we  differentiate  both  sides  of  equation  (C.l) 
with  respect  to  Xp,  the  X  coordinate  of  the  p   nucleus. 
Then  if  we  forget  that  the  spin-orbitals  u . ,  and  hence  the. 
charge  density  p,  depend  on  the  nuclear  positions  implicitly 
we  will  get : 


9<E 


xor   - 


*XP 


i-    f  p(?  )  _1_  (      2Zp      )  dV  -t-i-(    2ZP 


(C.3) 

which  is  what  the  Hellmann-Feynman  Theorem  contains  if  u-^ 

and  <E>  are  derived  from  the  exact  hamiltonian  instead  of 

the  Xa  hamiltonian  (C.l).  We  may  wonder  that  the  fact  that 

we  forgot  the  dependence  of  xi^   on  Xp  would  cause  some  errors. 

However,  we  show  in  the  following  that,  in  the  case  that  a 

is  constant  everywhere,  the  terms  arising  from  this  omission 

just  cancel  each  other,  and  thus  equation  (C.3)  exactly  holds. 

To  show  this,  we  consider  the  omitted  terms  F   .,: 

omit 


F 
omi 


+  I  il^T  ^ll  *\   ♦  |  \   «l3P1p(f1)  ]d3r2^ 


V 
2   9P(r2) 

9XP 


»J 


i3-  „(;  i  t  ri3.   2    3p(ri) 

2P(r2)  j  d  r1-?l2^-- 
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I  «(3/401/3  J  [^-(P(?)4/3)d3ri 


Since    p(r.)  =  2-  n.u*(r)u.(r) 
j   J  J     3 

therefore    !fi£l  -  T  (  ^-  ^   1°  ll£  ) 
3xp      .    3ui  3XP   3U*  3Xp 

■  -  X  <  niul(?)|^  ♦  niUi(r,^i  )   (C.4) 
=   T  n  (  3ui'?l>  f   v2  +  f-    ~22P 


Hence   F 

om 


12 


9         1/3  ■*.   1/3      -v     3 

-  I  cx(3/4tt)    (4/3)ps  (r]L)    }  u.(rx)  d^ 


+   (complex  conjugate) 


We  note  that  the  term  inside  the  bracket  together  with  u  (r  ) 

i   1 

in  the  above  equation  is  just  as  the  following: 

1       1 


2    <r-    -22p      r 


2Zp    _   r  2p(r2)  a  3    +  7 


d°r_  -  6o(3/4tt)  pQ  (r- )   } 


£?   l^-S   I  '   i    r_  **  *2  "  ^ '-^s,-l 
P-1   '  1   p  ' 


12  1 


x  ui<r1) 


From  equation  (2.4),  we  see  that  it  is  just  equal  to  £.u.(r_) 

ill 


Hence,  we  have: 
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F     =  V  n-  I  l^i^il  e.u.(r.  )d3r-  +  (complex  conjugate) 
omit     •     J    3X     li-L    -l 

=  Vne    3    f  u*(?1)u.(?1)d3r 
Z->   i  i  gXp   J  i   1   l   1    1 

On  account  of  the  normalization  of  u^,    we  have: 


_3 
3X 


{  [  uj(r1)ui(r1)ddr1  }  =  0 
P    J 


Hence,    FQmit  =  0 

One  small  point  about  this  derivation  is  that  in 
equation  (C.4),  we  had  implicitly  assumed  the  value  a  is 
constant  everywhere  so  that  ~-  =   0  and  so  —  does  not 
involve  the  term  |g  .  This  assumption  of  constant  a  does 
not  cause  any  trouble  if  we  are  considering  a  molecule  formed 
from  a  single  type  of  atoms,  such  as  a  homonuclear  diatomic 
molecule.  However,  some  corrections  have  to  be  made  in  the 
treatment  of  a  molecule  formed  from  different  types  of  atoms. 


APPENDIX  D 


PROOF  OF  THE  INDEPENDENCE  OF  DIPOLE  MOMENT 
ON  THE  CHOICE  OF  THE  ORIGIN  FOR  A  NEUTRAL  MOLECULE 


In  Chapter  III  and  Chapter  V,  we  claim  that  for  a 
neutral  molecular  system  the  value  of  the  dipole  moment  is 
independent  of  the  choice  of  the  origin  of  the  coordinates, 
so  that  we  can  locate  the  origin  at  the  center  of  the  outer 
sphere.  We  shall  prove  it  in  the  following: 

The  dipole  moment  y  computed  with  reference  to  an 
origin  0  is  defined  as: 


y >  X  ZA  ~   J  p(?)?  d3r  (D-1} 


where   p(r)  is  the  electronic  density  distribution, 
Z'-  are  the  atomic  numbers  of  the  nuclei, 

R.:  are  the  position  vectors  of  the  nuclei  with  the 

■  -*■ 
origin  at  0,  similarly  for  the  electronic  coordinate  r. 

Now,  suppose  the  origin  of  the  coordinate  system 

is  removed  to  0 '  by  the  transformation  T  (figure  D.l).  Then 

the  dipole  moment  y'  computed  with  reference  to  0'  is: 


y'  = 


^T  ZiR!  -   f  p(r')r'  d3r'  (D.2) 

i  * 
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Fig.  D.l  Transformation  of  Coordinates 


14 1 


where  R!  and  r'  are  measured  with  respect  to  the  origin  0 
From  figure  (D.l)  we  see: 


R^  =  Ri  -  T  and       r'  =  r  -  T 


Substituting  in  equation  (D.2),  we  have: 
V  =   2-  Zi(Ri-T)  -   J  p(r')(r-T-)dV 

=  {  ZiRi  -  j  p(r')r  d3r'}  -  {  £  Z.  -   j  p(r')d3r'}  T" 

-*■  -*■ 

Since  r'  and  r  are  referring  to  the  same  point,  we 

->■     ->• 
have  p(r')=p(r).  Hence  the  term  inside  the  first  bracket  is 

-*■ 
equal  to  p,  the  dipole  moment  value  computed  with  respect  to 

0.  Also,  the  term  inside  the  second  bracket  is  equal  to  the 

net  charge  of  the  system.  Since  we  have  a  neutral  system, 

the  net  charge  is  equal  to  zero.  Thus,  we  have: 

y'  =  u 

Since  T  is  arbitrary,  we  have  proved  that  the  dipole  moment 
is  independent  of  the  choice  of  the  origin  for  a  neutral 
system. 


APPENDIX  E 


EXPANSION  THEOREMS  FOR  THE  BESSEL  FUNCTIONS 
AND  THE  NEUMANN  FUNCTIONS 


Starting  from  the  well  known  plane-wave  expansion: 


e     =  4  tt   2,,  x   J0(kr)  Y0(r)  Y0(k> 


where  we  use  real  spherical  harmonics.  Then  we  express  the 
following  plane  wave  in  the  two  equivalent  expansions: 

eik.(r2-V  =  4,  £  irj£,(k|?2-?il)Y^(?2-?l)Yr^)  (E.l) 

i£. (?--?.)    iS.?2   -i£.r\. 
e     21=e^e     x 

V"   £  '         m1  -»-   m'  -► 
=-4tt  2_^ix  j£,(kr2)YJi.(r2)Y£l(k) 


i  > 


•£" 


4irT  i    j   (kr  )Ym  (r  )Ym  (S)       (E.2) 

'  V.    r   l  £'*  l  2." 

Multiplying  both  sides  of  equation  (D.l)  by  Y  (k)  and  integ- 
rating over  the  solid  angle  Q,,    we  have: 


eiS-(?2-?l>Y™(J)dnk  =  ^iV^V^l'^VV      (E.3) 
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by  the  orthonormallty  of  the  spherical  harmonics.  Similarly, 
for  (E.2) ,  we  have: 

J  eiS-<V?l>Y^(S)dnk 

=  (4TT)"5  X   X  i     I(£'m';£"m";£m)j  ■,(kr2)Y™,  (r?) 

x  J^.Ckr^Y^C^)  (E.4) 


Comparing  (E. 3)  with  (E.4),  we  find: 

j£(*iV?ii>^<v*i) 

=  ^  X  X  ^  ~£"~£  I(^m';£"m";ilm)  j  ,  (kr  )Y™',  (r  „) 

m"  ->- 
j   (kr  )Y   (r  )  (E.5) 

0"    1   5,'    1 


m  .-»■ 
which  is  the  expansion  for  j  (kr)Y  (r). 

To  derive  a  similar  expansion  for  n  (kr)Y"1(r),  we 

use  the  well  known  spherical -wave  expansion: 


e^'^l'/lr^-rj  =  4*  ik  ^  ^(kr^j^kr^Y^?^^?^ 

for   r1>  r2       (E.6) 
where  h*(x)  is  the  spherical  Hankel  function  of  the  first 
kind  and  is  related  to  j„  and  n«  by  h^(x)=jji(x)  +  inn  (x) . 
Thus ,  we  have : 
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ik1?n-?,-?. 

2   *1   x3 


e  *      -1      J   /  \r9-r~  -r 


4Trik  Th^(k|r2-r1|)j£(kr3)Y^(?2-?1)Y1J(?3),     |?  ^?    |  >r      (E.7) 
4,1k    ShJCkr^J    (kl?^!)^?^^?^),    r^lV^I    (E'8) 


->-  -*■   v,Tin,  -*•  •*■ 


With  the  expansion  for  j  (k|rQ-r„|)Y  (rQ-r  )  as  expressed 
(E.5),  it  follows  from  (E.8)  that 


eik|r2-r1-r3|/|?2_?i_?3 


-4irik  X  ^(kr  )Ym(r"  )4tt  T  £  1*'"*""* 


I(£'m'  ;£"m";lra) 


,m' 


■  x  j   (kr  )Y   (r  )j   (kr  )Y   (r  )      if  r  >|r  -r  I     (E.9) 
J£'    2  I'      2      1"         3  H"      3  1  '  2   31 

If  r- ,  r„,  r_  satisfy  both  conditions  in  (E.7)  and  (E.9), 
i .e. I r_-r1 I >r„  and  r1>|r„-r„|  then  we  can  combine  (E.7)  and 
(E.9)  together,  and  after  rearranging  the  dummy  summation 
index  in  (E.9)  and  equating  the  coefficients  of  ,ip  (kro)Y.  (rJ 
we  have  the  following  identity: 

bJCklr^OYjc?^) 

hJl„(kr1)Y^t,(r1)  ^  i        I(£  'm  '  ;  Zm;  l"ml  ) 

i'vn"  JJ'V 

x  J£-(kr2)Y£'(r2)     if  ,r2"ril>r3  and  ri>|r2~r3'     (E'10) 
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Since  h£,=Jp+ino  >  and  Jo  satisfies  the  same  expansion  theorem 
as  in  (E.10),  thus  we  have  the  expansion  theorem  for  the  Neu- 
mann function: 


n  <k|r  -r"  |)Ym(?  -r"  ) 
£     2   1    Z      2   1 

=  4tt  T  ^i*'"*"*"  i(jL.m.;£m;£"m")n£n(kr1)Y™"(r1) 

x  Ji,(kr2)Ym,(r2)    if  |r2-r1|>r3  and  r1>|r2-r3|     (E.ll) 

->      -*■  -»- 

Since  (E.ll)  involves  only  r-^  and  r2,  irrespect  of  rg,  thus 

-»■ 
we  can  take  r,  to  be  located  at  the  origin,  and  then  the  con- 
ditions for  (E.ll)  reduce  to   r1>r2. 

For  the  case  r..<r2,  we  can  derive  the  expansion  for 
the  Neumann  functions  in  a  similar  way  by  rearranging  the 
arguments  in  (E.7)  and  (E.8).  The  result  is: 

=  47r  £  X„i^~£~£"  W^' ;im;£"m")  JJl„(kr1)YJi„(r1) 

m'  -*■ 
x  nA-,Ck"r2)Y£,(r2)        if   r   <  rg  (E.12) 


APPENDIX  F 


EXPANSION  THEOREM  FOR  Y?(r")/r2 

1   A    A 


With  the  expansion  for  the  Neumann  function  (Appendix 
E ) ,  we  have : 

n  (k|x.-R.  ■.  |  )Y™(x.-R.  .) 

=  4tt  V   X  j/'-*-*"  l(£"m";^mn£,m')  n   (kR   )Ym"(S   ) 
i>'  *>"  a"    ^   r'   ij 

m'  -*■■  ■+  -*■ 

x  j  ,(kx.)Y  ,(x.)        if   x.  <  R. .  (F.l) 

Since,  as  k'  +  0,  we  have: 


fc+1 
D  (kx)  +  -  (2fc-l)! !/(kx) 


and  Jo(kx)  "*  "  (kx)£/(2£+l)!  ! 


where   n!!  =  n(n-2) (n-4)  . .  . 1 

and    (-1)! !  =  1 

Thus,  as  k  ■■*■   0,  equation  (F.l)  becomes 
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-■ { ' (2£-l) ! !/(k|x.-S. .| >A+1  >  Y^x.-R..) 
i   ij         £   l   ij 


i1'-1-*"  IC£"m";£m;£-m-)  {-  ^'^V'   Y^R..) 

(kxn-  )     mi  -»- 

x  — 3J Y„,(x,  )}  if   x.  <  R.  , 

(2£'+l)!  !   A*   i'J  i     1J 

£+1 
Multiplying  k    on  both  sides  and  taking  k->-0,  then  the  only 

nonvanishing  term  on  the  RHS  is  the  one  with  £"=£'+£,  and  we 

have: 


m  -*■     ■+      v  ,  .-»■  -»•  ■  .  £+1 
Y  (x.-R.  .)/  x.-R.  . 
A   i  .  ij     i   ij 

-'41.  1  Z  i-"i(fw,.";i,.H'.»')  »    121^2izllil 

«C*   W  (2*-l)!!(2l'+!)!! 

xf/Rj:+l+1y™;(;.)Y™';+J(R..)     ifx.<R..  (F.2) 


Now,  in  the  treatment  of  diatomic  molecules,  we  have  centers 
A  and  B.'  Thus,  let 


R .  .  —  R   —  R_ 
ij     A     B 


x.  =  r   +  R.  .  =  rn 
l  ■    A     ij     B 


then   x.  -  R. .  =  r. 
l     ij     A 


Substituting  in  (F.2)  and  considering  the  case  for  £=1,  m=0J 
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Yl(rA)/i\  =  4tt  V   A  i"2  l(£'+l,m";l,0;£\m')  (2£'+1)!  ! 
A   A      /w   *M  (2£'+l)M 


*\ 


£ '  ,„£  '+2  „m'  ,-»-  x„m" 


*  rB  /RAB   V'W+l'W     lf  rB<RAB 


Since  RA-RB  is  along  the  z-axis,  therefore 


m"        ■*■     ■*■  0    -*-  ->■ 

Y  ,   (R  -R  )  =  6     Y     (R  -R  ) 
l'+l      A   B     m"0   r+1   A   B 


which  implies  that  the  only  nonvanishing  term  in  the  summa- 
tion over  m'  is  the  one  with  m'=0.  Thus,  we  have: 


C0 

Y  (r  )/r   =  -  4tt  >   I(  I '  +1 ,  0;  1 ,  0;  I '  ,  0)Y    (R  -R  )/R 
1   A    A         fc*„  fc'+l   A   B  '  AB 


+2 


t'~Q  *•      T-L    A    D      AB 

'b'-YJ-(V-        "  VRAB  (F-3) 


Similarly  we  can  prove  the  expansion  theorem  for  the  case 
r   >  R   by  using  (F.12).  The  result  is: 

Y°(?A)/r  -  4,  V  Kl'-l.Oil.Ojf.OjR^-HJ,   (i  -5  ) 

*Y!.(?0)/ro'+1      if  VRA0  <F'4> 


APPENDIX  G 

EVALUATION  OF  THE  SINE  AND  COSINE  INTEGRALS 

In  the  force  and  dipole  moment  calculations,  we  have 
to  evaluate  the  following  integrals: 

_J1H_   dt 

Jti  tP(a-t)^ 

where     f(t)  =  1,  cost,  sint,  and  p,q  are  integers. 

It  can  be  proved  that  the  evaluations  of  these  kinds 
of  integrals  wind  up  with  the  evaluations  of  the  sine  and 
cosine  integrals,  namely: 

J  cosx/x  dx  and  sinx/x  dx 

as  we  shall  show  below: 

(1)  q  <  0  • 

In  this  case,  we  can  use  the  binomial  expansion  for 

-q 
(a-.t)      •    Then,    we  have: 


-9    ,    i-,k      -Q-k  -+-<1-P 


Tl  sq  ^k  a  - "  * 


tp(a-t)q  fci      k 

Thus,      f    f(t)    =  y  c:q(-i)ka-q-k  ftA  mit 

Jti   tP(a-t)^         t^o     k  Jt.    tp+q 


1*9 
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(2)  p  <  0 

In  this  case,  we  change  the  variable  t  to  x:  x=a-t 
Then,  we  have: 

cost  =  cos(a-x)  =  cosa  cosx  +  sina  sinx    etc. 


and  tP(a-t)q  =    (a-x)PxQ 


r   cost    dt  , 

JtitP(a-t)^  Jfc.ti 


A-ta 

cosa 

xH(a-x)J 


Thus     I    —   dt  =   \     cosa  cosx  +  sina  sinx  dx  etc 


which  can  be  evaluated  just  like  Case  (1). 
(3)  q  >  0,  p  >  0 

Integrating  by  parts  and  using  the  following  identity 

D+a-z  I        Q_+     D-S-l 

dt 


f  dt  l-p-qr       V  p+q-2  f     a-t ,p-s-l 

J  ?w?  =  a      {  h  Cs     /(p-s_1)  J  (") 

>i*p-i 

Cp-1        logl"T-|    > 


we  wind  up  with  the  integrals  encountered  in  Case  (1)  or 
Case  (2)  according  to  whether  p-s-1  is  positive  or  negative. 
Thus,  for  all  cases,  we  winds  up  with  the  evaluation 
of  integrals : 


f*1      m  f** 

cosx/x  dx 

"X,  J%. 


sinx/x  dx 
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For  negative  values  of  m,  these  can  be  reduced  to 
the  evaluation  of  sine  and  cosine  functions  by  the  formula: 

f  -m     ,     -m  .    ^    f   -m-1     , 
I  x  cosx  dx  =  x  smx  +  m  J  x    sinx  dx 

f  -m  .  -m  f  -m-1 

J  x  sinx  dx  =  -x  cosx  -  m   I  x    cosx  dx 

For  positive  values  of  m,  we  have  similar  relations 
j  cosx/xm  dx  =  -  cosx/  (rn-Dx"1-1  -  ^  j  sinx/x^-l  dx 

f  .   /  ni  ,  ,  ,      -  ,  m-1     l   f     ,  m— 1  , 

j  sinx/x  dx  =  -  sinx/  (m-l)x     +  — —  J  cosx/x    dx 

(m  ^  2  ) 

which  reduces  to  the  evaluation  of  the  integrals: 


cos/x  dx  and  I  sinx/x  dx 


P 


*, 


These  two  integrals  can  be  evaluated  as: 

I  cos/x  dx  =  ci(x  )  -  ci(x  ) 
ht  2  X 

1  sin/x  dx  =  si(x2)  -  si(x1) 
*. 

where   ci(x)  and  si(x)  are  defined  by:  . 
ci(x)  =  C  +  log  x  +    I  (cost-l)/t  dt 


152 


si(x)  =    J  sint/t  dt  -  tr/2 


where  C  is  the  Euler   constant  0.57..... 

The  functions  ci(x)  and  si(x)  can  be  computed  with  a 
subroutine  called  SICI  which  is  stored  in  the  IBM  Scientific 

Subroutine  Package: 

2 

For  x£4  and  t=l6-x  : 

{si(x)+7T/2}/x  =  4.395509E-1  +  1.964882E-2  t  +  6.939889E-4  t2 
+  1.374168E-5  t3  +  1.568988E-7  t4 
+  1.753141E-9  t   +  e 
with  |e|  <  8.2E-8 
{-ci(x)-C-log|x| }/x2  =  1.315308E-1  +  4.990920E-3  t 

+  1.185999E-4  t2  +  1.725752E-6  t3 
+  1.584996E-8  t4  +  1.386985E-10  f5  +  e 
with  |e|<  5.6E-9 
For  x>4  and  t=4/x: 
si(x)  -  -  (4/x)  {cosx  P(x)  +  sinx  Q(x)   +  e} 
ci(x)  =  (4/x)  {sinx  P(x)  -  cosx  Q(x)   +  e} 
where  | e |  <  2.3E-8  x~ 

and  P(x)  =  2.500000E-1  -  6.646441E-1  t  -  3.122418E-2  t 

3  4  5 

-  3.764000E-4  t   +  2.601293E-2  t   -7.945556E-3  t 

-  4.400416E-2  t6  +  7.902034E-2  t?  -  6.537283E-2  t** 

9  10 

+  2.819179E-2  t   -  5.108699E-3  t 

Q(x)  =  2.583989E-1  +  6.250011E-2  t  -  1.134958E-5  t2 

-  2.314617E-2  t3  -3.332519E-3  t4  +4.987716E-2  t5 

6  7  8 

-7.261642E-2  t   +  5.515070E-2  t   -2.279143E-2  t 

+  4. 048069E-3  t9 
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